{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import networkx as netx\n",
    "import numpy as np\n",
    "import DaPy as dp\n",
    "from random import randint, random"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " - read() in 0.037s.\n"
     ]
    }
   ],
   "source": [
    "data = dp.read('csv/Updates_NC.csv').data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def decode(val, code='utf-8'):\n",
    "    if isinstance(val, bytes):\n",
    "        return val.decode(code)\n",
    "    return val\n",
    "data = data.map(decode, cols=['报道时间', '省份', '城市', '消息来源'], inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " 报道时间 | 新增确诊 | 新增出院 | 新增死亡 \n",
      "----------+----------+----------+----------\n",
      " 1月17日  |    17    |    4     |    0     \n",
      " 1月18日  |    4     |    3     |    0     \n",
      " 1月19日  |    17    |    4     |    0     \n",
      " 1月20日  |    77    |    5     |    1     \n",
      " 1月21日  |    60    |    0     |    2     \n",
      " 1月22日  |   105    |    3     |    8     \n",
      " 1月23日  |    62    |    0     |    17    \n",
      " 1月24日  |    70    |    3     |    6     \n",
      " 1月25日  |    77    |    1     |    15    \n",
      " 1月26日  |    46    |    8     |    7     \n",
      " 1月27日  |    80    |    2     |    18    \n",
      " 1月28日  |   892    |    3     |    22    \n",
      " 1月29日  |   315    |    26    |    19    \n",
      " 1月30日  |   356    |    7     |    25    \n",
      " 1月31日  |   378    |    21    |    30    \n",
      "  2月1日  |   576    |    36    |    33    \n",
      "  2月2日  |   894    |    32    |    32    \n",
      "  2月3日  |   1033   |    53    |    41    \n",
      "  2月4日  |   1242   |    79    |    48    \n",
      "  2月5日  |   1967   |    65    |    49    \n",
      "  2月6日  |   1766   |    58    |    52    \n",
      "  2月7日  |   1501   |   103    |    64    \n",
      "  2月8日  |   1985   |   164    |    67    \n",
      "  2月9日  |   1379   |   179    |    63    \n",
      " 2月10日  |   1920   |   167    |    73    \n",
      " 2月11日  |   1552   |   162    |    67    \n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "wuhan = data.query(' 城市 == \"武汉市\" and 报道时间 > \"1月16日\"').reverse()\n",
    "wuhan = wuhan.groupby('报道时间', max, apply_col=['新增确诊', '新增出院', '新增死亡'])\n",
    "wuhan.insert_row(0, ['1月17日', 17, 4, 0]) # 补充两条国家卫检委数据\n",
    "wuhan.insert_row(6, ['1月23日', 62, 0, 17])# 补充两条国家卫检委数据\n",
    "wuhan.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### SEID模型\n",
    "易感数量 Susceptible：可能感染该病的健康人员数量\n",
    "\n",
    "潜伏数量 Exposed：已感染病毒，但仍未表现出症状\n",
    "\n",
    "患者数量 Infected： 已感染病毒，并表现出了症状\n",
    "\n",
    "死者数量 Dead：死亡患者人数\n",
    "\n",
    "治愈数量 Cured：治愈患者人数\n",
    "\n",
    "潜伏感染率 $\\alpha_{1}$：由潜伏者导致的传染概率\n",
    "\n",
    "患者感染率 $\\alpha_{2}$：由患者导致的传染概率\n",
    "\n",
    "确诊率 $\\beta$：潜伏者表现出症状的概率\n",
    "\n",
    "治愈率 $\\sigma$：患者治愈的概率\n",
    "\n",
    "死亡率 $\\gamma$：患者死亡的概率\n",
    "\n",
    "\n",
    "$$\n",
    "\\frac{dS}{dt} = - \\alpha_{1} E - \\alpha_{2} I + \\sigma I \\\\\n",
    "\\frac{dE}{dt} = \\alpha_{1} E + \\alpha_{2} I - \\beta E\\\\\n",
    "\\frac{dI}{dt} = \\beta E - \\sigma I - \\gamma I\\\\\n",
    "\\frac{dD}{dt} = \\gamma E\\\\\n",
    "\\frac{dC}{dt} = \\sigma I \\\\\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import dual_annealing, minimize\n",
    "from sklearn.metrics import r2_score\n",
    "from collections import namedtuple\n",
    "from matplotlib import pyplot as plt\n",
    "plt.rcParams['font.sans-serif'] = ['SimHei']\n",
    "\n",
    "SEIR_PARAM = namedtuple('SEIRparm', ['alpha1', 'alpha2', 'beta', 'sigma', 'gamma'])\n",
    "\n",
    "class SEIR(object):\n",
    "    def __init__(self, P=None):\n",
    "        self.P = P\n",
    "        \n",
    "    def _forward(self, S, E, I, D, C, param, max_iter):\n",
    "        a1, a2, b, s, g = param\n",
    "        est = dp.Table(columns=['S', 'E', 'I', 'D', 'C'])\n",
    "        for t in range(max_iter):\n",
    "            S_ = S - a1 * E - a2 * I + s * I\n",
    "            E_ = E + a1 * E + a2 * I - b * E\n",
    "            I_ = I + b * E - s * I - g * I\n",
    "            D_ = D + g * I\n",
    "            C_ = C + s * I\n",
    "            S, E, I, D, C = S_, E_, I_, D_, C_\n",
    "            est.append_row([S, E, I, D, C])\n",
    "        return est\n",
    "    \n",
    "    def _loss(self, obs, est):\n",
    "        assert len(obs) == len(est)\n",
    "        loss = ((np.log2(obs + 1) - np.log2(est + 1)) ** 2).sum()\n",
    "        self.lossing.append(loss)\n",
    "        return loss\n",
    "    \n",
    "    def _optimize(self, param, s, e, i, d, c, obs):\n",
    "        est = self._forward(s, e, i, d, c, param, len(obs))\n",
    "        return self._loss(obs, est['I', 'D', 'C'].toarray())\n",
    "    \n",
    "    def fit(self, initS, initE, initI, initD, initC, Y):\n",
    "        self.lossing = []\n",
    "        args = (initS, initE, initI, initD, initC, Y['确诊', '死亡', '治愈'].toarray())\n",
    "        param = [(0, 1),] * 5\n",
    "        # result = minimize(self._optimize, [random()] * 5, args=args, bounds=param)['x']\n",
    "        result = dual_annealing(self._optimize, param, args=args, seed=30, maxiter=10)['x']\n",
    "        self.P = SEIR_PARAM(*result)\n",
    "    \n",
    "    def score(self, initS, initE, initI, initD, initC, Y, plot=False):\n",
    "        est = self.predict(initS, initE, initI, initD, initC, len(Y))['I', 'D', 'C']\n",
    "        loss = self._loss(Y['确诊', '死亡', '治愈'].toarray(), est.toarray())\n",
    "        est.columns = ['确诊', '死亡', '治愈']\n",
    "        r1 = r2_score(Y['治愈'], est['治愈'])\n",
    "        r2 = r2_score(Y['死亡'], est['死亡'])\n",
    "        r3 = r2_score(Y['确诊'], est['确诊'])\n",
    "        if plot:\n",
    "            self.plot_predict(Y, est)\n",
    "            print(' - 平均潜伏期为：%.2f天' % (1.0 / self.P.beta))\n",
    "            print(' - 病毒再生基数：%.2f' % (self.P.alpha1 / self.P.beta + (self.P.alpha2 / self.P.sigma + self.P.alpha2 / self.P.gamma)/ 2))\n",
    "            print(' - 确诊R2：%.4f' % r3)\n",
    "            print(' - 死亡R2：%.4f' % r2)\n",
    "            print(' - 治愈R2：%.4f' % r1)\n",
    "            print(' - 模型R2：%.4f' % ((r1 + r2 + r3) / 3))\n",
    "            print(' - 模型总误差：%.4f' % loss)\n",
    "        return loss, (r1 + r2 + r3) / 3\n",
    "    \n",
    "    def plot_error(self):\n",
    "        plt.plot(self.lossing, label=u'正确率')\n",
    "        plt.legend()\n",
    "        plt.show()\n",
    "    \n",
    "    def plot_predict(self, obs, est):\n",
    "        for label, color in zip(obs.keys(), ['red', 'black', 'green']):\n",
    "            plt.plot(obs[label], color=color)\n",
    "            plt.plot(est[label], color=color, alpha=0.6, label=label)\n",
    "            plt.legend()\n",
    "            plt.show()\n",
    "            \n",
    "    def predict(self, initS, initE, initI, initD,initC, T):\n",
    "        return self._forward(initS, initE, initI, initD, initC, self.P, T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "c:\\python36\\lib\\site-packages\\ipykernel_launcher.py:28: RuntimeWarning: invalid value encountered in log2\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "潜在患者：25.0000 | R2：0.8450 | 误差： 90.591710\n",
      "潜在患者：100.0000 | R2：0.8450 | 误差： 90.591710\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD6CAYAAACoCZCsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XmYVNWZx/HvK3Szi2wii4g6iBgVgdaAaIILCGqMOCoaRuMy0SwmTjSJOubJjMZlQhJn1KAJCa5xTNBExVEUkEUEBLtZFWIQQWwXaGnW7gZ6OfPHW51qoJfqrW511e/zPPep6lO3us6lqfOee1YLISAiIpnnkKgzICIi0VAAEBHJUAoAIiIZSgFARCRDKQCIiGQoBQARkQylACAikqEUAEREMpQCgIhIhmoddQZq071799C/f/+osyEi0qLk5eV9EULoUdd5KR0A+vfvT25ubtTZEBFpUczso0TOUxOQiEiGUgAQEclQCgAiIhkqpfsAqlNaWkp+fj579uyJOitNrm3btvTt25esrKyosyIiGaDFBYD8/Hw6depE//79MbOos9NkQghs3bqV/Px8jj766KizIyIZoMU1Ae3Zs4du3bqlVeEPYGZ069YtLe9sRCQ1tbgAAKRd4V8pXa9LRFJTiwwAIiJp7Y03YPnyZv8YBYAm8tJLL7F06dJaz/niiy/Ytm1bknIkIi1SeTlMnw6rVzf7R9UZAMyss5nNMLOZZvaCmWWb2SYzmxc7Toqdd5eZvWNmk6u8N6G0lubTTz/l6quv5pprruGqq65i48aN3HTTTUyZMoUf/ehHfPzxxwBcfvnl5OXl/eN99913H+PGjYsq2yLSEqxbB3v2wMknN/tHJXIHMBF4IIQwBvgcuB14NoQwKnasNrNhwBnAacAWMzs30bTmuKjm1qtXLx566CG2bdvGE088wY9+9CMmT57MQw89xJo1a+jbty8A999/P7/73e8A+Oyzz1i2bBmjR49m2rRpUWZfRFLZqlXQujUMGtTsH1VnAAghPBJCmBX7sQdQBlxoZkvNbKqZtQa+CvwlhBCA14Ez65HW4pgZixcvZvTo0ZSVlXHqqafywAMPMGLECNauXctZZ53FkiVLOPbYY5kyZQolJSVMnDiRSZMmMWHCBB5//HFuueUWNm7cGPWliEgqCcEDwMCB0KZNs39cwvMAzGwE0AWYBTweQvjMzJ4Czgc6AOtjpxYCPfFAkUjagZ9zA3ADQL9+/WrP1LRpEGtuaTJHHgmXX17naVOnTiWEwNy5c5k0aRIlJSUcccQRHH/88eTm5jJ79mxuvvlmhgwZwqZNm7j22mtZuHAheXl5PPfcc0yaNIlDDlEXjIhUsXkzFBTAuclpHEkoAJhZV+Bh4J+Bz0MIe2Mv5QIDgN1Au1haR/zOItG0/YQQpgBTAHJyckL9Lic5Zs+ezfTp01mzZg2333475eXlPPHEE2RnZ9O+fXt27drF1KlTGTlyJDNmzOCFF17g3//93yksLORXv/oVZ599NvPmzaN9+/ZRX4qIpJJVq/wxCe3/kEAAMLNs4DngjhDCR2Y2zczuBd4FLgbuA/YBlwN/AgYDG4G8BNMaLoGaenPo168f11577T9+Li4uZuzYsXTv3p3+/fuzZs2af7xmZmRnZzNu3DiOOuooxo8fz913382cOXMYM2YM2dnZUVyCiKSiVaugb1/o2jUpH5dIG8T1wFDgTjObB7wHPA2sABaHEGYDbwFDzOxBYp3E9UhrcY477ji6dOnCxIkTyc3NZcOGDbRv357i4mJ2795NSUnJP87Ny8vj/vvvZ8OGDVx22WVMnjyZ/Px87rnnHg0JFZG4oiJYvz5ptX9I4A4ghPAo8OgByXcdcE5FbETPBcCDIYQNAImmtUQVFRU888wzPProo7z00kvcdtttzJ8/n+OPP57S0lIACgsLWblyJb///e857LDDOOecc7j33nvZvn07c+bMUROQiMS9+y5UVKRWAEhUCKEEeL4haS3ReeedR/fu3bntttt47LHHGDRoEGVlZWRlZZGXl8eECRMoLS3lueeeo0+fPowfP5527dpx4403cm6SOnhEpAVZvRo6dYIkboNrPiIzNeXk5IQDt4Rcu3Ytg5IwPrap7d27lzYJDOtqqdcnIo1QXg633gpDhsA3v9noX2dmeSGEnLrOa5HjEFM5aNUkkcK/JV6XiDSBDz6AkpKkNv9ACwwAbdu2ZevWrWlXWFbuB9C2bduosyIiyVY5+/eEE5L6sS1uQ5i+ffuSn59PQUFB1FlpcpU7golIhkni7N+qWlwAyMrK0o5ZIpI+Nm+GLVvgnHOS/tEtrglIRCStrFzpjyedlPSPVgAQEYnS6tU++7dbt6R/tAKAiEhUiot9BFAEtX9QABARiU4Es3+rUgAQEYnKqlU++zeigS0KACIiUSgvh/fe8+Yfs0iyoAAgIhKF9eu9DyCi5h9QABARicbKlZHM/q1KAUBEJAqrV0cy+7cqBQARkWTbvNmPiIZ/VlIAEBFJtiTv/VsTBQARkWRbtQr69Ilk9m9VCgAiIslUOfs34to/KACIiCTXe+9FOvu3KgUAEZFkqpz9m8S9f2uiACAikiwVFb7+z4knwiHRF7/R50BEJFN88EHks3+rUgAQEUmW1asjn/1blQKAiEiyrFwJxx0HbdtGnRNAAUBEJDm2bPHZvynS/AMKACIiyZEis3+rUgAQEUmGVaugd+/IZ/9WpQAgItLcioth3bqUqv2DAoCISPNLodm/VSkAiIg0t9WroWPHyPb+rUmdAcDMOpvZDDObaWYvmFm2mU01s8Vm9tMq5zU4TUQkbZWX++zfk05Kidm/VSWSm4nAAyGEMcDnwBVAqxDCCOAYMxtgZpc0NK15LktEJEW8+y4UFcHQoVHn5CCt6zohhPBIlR97AP8C/E/s55nAGcAQYFoD09Y1PPsiIilu0SI49FBf/yfFJHw/YmYjgC7Ax8AnseRCoCfQoRFpB37ODWaWa2a5BQUF9boYEZGUsmuXD/8cPjzlmn8gwQBgZl2Bh4HrgN1Au9hLHWO/ozFp+wkhTAkh5IQQcnr06FHf6xERSR1Ll/ronxEjos5JtRLpBM4GngPuCCF8BOThTTcAg4GNjUwTEUk/IXjzT//+PgEsBdXZBwBcDwwF7jSzO4HHgavMrDcwDhgOBGBBA9NERNJPfr4fV14ZdU5qVOcdQAjh0RBClxDCqNjxJDAKeBs4K4SwI4Sws6FpzXFRIiKRW7jQl34+7bSoc1KjRO4ADhJC2EZ8NE+j00RE0kpZmbf/n3IKtG8fdW5qlHrd0iIiLd2qVT72//TTo85JrRQARESa2uLFcNhhMGhQ1DmplQKAiEhT2rnTZ/+m6Nj/qlI7dyIiLc3bb/vY/xRv/gEFABGRplM59v+YY6DnQQsdpBwFABGRpvLRR/DZZ42v/S9bBp98Uvd5jaQAICLSVBYtgqwsyMlp+O9YsgTOPhuuv77p8lUDBQARkaZQWgrvvANDhkC7dnWfX53Fi2HMGN83eMqUps1fNRQARESawsqVvvdvQ5t/Fi2C886Dww+H+fOhX7+mzV81FABERJrCokXQpQsMHFj/9771lhf+RxwB8+ZB375Nnr3qKACIiDTW9u2wZo0v+1zfsf9vvgljx0KfPl749+nTLFmsjgKAiEhjvf22DwGt77r/8+bBuHFw5JEwd27Sl41WABARaYzKsf8DBnj7faLeeAPOP9/3C5g3D3r1aq4c1kgBQESkMTZsgM2b61f7nzULLrwQjj3Wa/4RTRpTABARaYxFiyA7G4YNS+z811+Hr30NjjsO5syp311DE1MAEBFpqH37fOz/0KHQtm3d58+YAV//uq8S+sYbEPG+5woAIiINtXw57NkDI0fWfe4rr8DFF8MJJ3jh37178+evDgoAIiINtXixz9odMKD2815+GcaPh5NOgtmzoWvX5OSvDgoAIiINUVgIf/ubz/w1q/m8xYvh0kt9e8hZs1Km8AcFABGRhqkc+z98eM3n5Od7zf/II+G113ymcApp0KbwIiIZrXLs/8CBNbfll5R4m39xsbf5p1DNv5LuAERE6uuDD6CgoOaF30Lw5ZyXLYNnnoEvfSm5+UuQ7gBEROpr0SJo08aXfq7OpEnw7LNw770+5j9F6Q5ARKQ+9u6FvDzf9KVNm4Nff+UVuOMOmDDBH1OYAoCISH0sXOhB4IwzDn5t7Vr4xjf8zuCxx2ofHZQCFABERBJVVgYzZ/q4/2OO2f+1bdvgoot8RvCLL0L79tHksR4UAEREErV0qRf0Y8fun15WBldc4ZvC//WvPuyzBVAnsIhIIioqfCz/kUcePKrnttv8zuAPf0hsWYgUoTsAEZFErFjhyz6PHbt/2/5TT8EDD8D3v+9DP1sQBQARkbqE4Ct5Hn64r/xZ6e234VvfgrPPhl//Orr8NVBCAcDMeprZgtjzPmaWb2bzYkePWPpUM1tsZj+t8r6E0kREUtqaNbBpk9f+K/f8/eQTX+ahb1+YNg2ysqLNYwPUGQDMrAvwJNAhlvRl4N4QwqjYUWBmlwCtQggjgGPMbECiac1zWSIiTahyHZ8vf9l/Linxwn/3bnjpJV8RtAVK5A6gHJgA7Iz9PBz4VzNbZmb3xdJGAdNiz2cCZ9QjTUQkda1fD3//O4weDa1be3PQDTf4RjB//COceGLUOWywOgNACGFnCGFHlaQZeEF+KjDCzE7G7w4+ib1eCPSsR9p+zOwGM8s1s9yCgoJ6X5CISJN67TXo0MEnfpWWekfvH/8IP/+57+7VgjWkE3hRCGFXCKEcWA4MAHYD7WKvd4z93kTT9hNCmBJCyAkh5PSIeLs0Eclw+fmwahWcc47P/r3gAnj8cfjZz+DOO6POXaM1JAC8bma9zKw9MAZ4F8gj3pwzGNhYjzQRkdT02mu+3s/AgfCVr/gm7lOnwl13pfwyD4loyESwu4C5wD7gtyGE983sM2CBmfUGxuH9BCHBNBGR1FNQALm5XviffbbPAH7lFTjvvKhz1mQSDgAhhFGxx7nA8Qe8ttPMRgGjgUmVfQaJpomIpJzXX4fPP/dlnTt0gDffrHn55xaqyZaCCCFsIz7Cp15pIiIpZft2ePJJmDvX7wBmzIB+/aLOVZPTTGARkapCgO9+F2bP9nH/b72VloU/aDE4EZG4sjK48UZv9jntNO/0rW7TlzShACAiAlBU5Es6/9//wSmnwF/+ktaFPygAiIj4Kp8XXuhbPZ59Nlxyia/xk+YUAEQks23Y4BO9Nm/28f35+TBuXNS5Sgp1AotI5ioq8uUctm2DWbN8cbfjjjt4u8c0pQAgIpmpclG3d9/1Tt+KCh/+ef75UecsaRQARCQzPfQQ/O//+qJuY8b4xK+jjoLjj6/7vWlCAUBEMs/8+XDrrd78c8cdsGwZbNnibf9psMZPohQARCSz5OfD5ZfDscf6fr5mPtP3iCN8+GcGUQAQkcyxdy9ceikUF8OLL8Khh/pyz/n5B2/2ngE0DFREMscPfgBLlvgkr0GDfIOXadO89n/aaVHnLul0ByAimeEPf4ApU+D2232iF8DMmfDFF3DlldCqVbT5i4ACgIikv6VL4Xvf831977nH07Zu9bb/YcMyauRPVQoAIpLetmyBf/5n6NXLx/tX1vSfe87b/C+7LNr8RUh9ACKSvsrKYMIEb+ZZuBC6dfP0996D5cvh4ouhS5do8xghBQARSV+33w7z5vnmLkOHelpZGfzpT3D44d4klMHUBCQi6enPf4Zf/xpuugmuvjqePmuWNwtdcQW0zuw6sAKAiKSf1avhuutg5EgPApUKC+HVV31v3y99Kbr8pQgFABFJL9u3w/jx0Lmzd/RmZ8dfe/55XwQugzt+q8rs+x8RSS87dsAFF8CmTd7236tX/LW1a33Dl69/Pd4ZnOEUAEQkPWzdCued50s7PPssnH56/LXKjt8ePTK+47cqNQGJSMv3+ecwapSv7f/iiz7uv6o5c/ycCRMgKyuSLKYi3QGISMv28ce+peOnn3oH79ln7//69u2+0fvJJ8NJJ0WTxxSlACAiLdf69V74b9vm6/pUbfap9PzzUF7utX/ZjwKAiLRMa9fCuefCnj3exDNs2MHnvP8+vPMOXHghdO+e/DymOAUAEWl5VqzwbRwPOcR39zrxxIPPKS/3jt9u3XytfzmIOoFFpGVZsgTOOgvatoUFC6ov/AHmzvV+AXX81kgBQERajvnzvdmna1cv/AcMqP68HTvg5Zc9OJx8cnLz2IIoAIhIy/D6675p+5FHeuF/1FE1n/uXv8RXAs2wbR7rI6EAYGY9zWxB7HmWmb1sZgvN7LrGpomI1OnFF+Gii2DgQL8L6N275nPXrfNmojFjfMVPqVGdAcDMugBPAh1iSd8H8kIII4FLzaxTI9NERGo2bZpv5D5kiI/26dGj5nOLi+Hxx73jd9y45OWxhUrkDqAcmADsjP08CpgWe/4mkNPItP2Y2Q1mlmtmuQUFBQlfiIikoVmzYOJEH98/a1btm7eEAE884RO/brhh/0XgpFp1BoAQws4Qwo4qSR2AT2LPC4GejUw78POmhBByQgg5PWqL9CKS3las8CUdBg2C6dOhUx0NBnPmwMqV/p7+/ZOSxZauIZ3Au4F2secdY7+jMWkiIvvbuNGbcDp39uUdDjus9vM3bPAZv6eccvBSEFKjhhTAecAZseeDgY2NTBMRiSss9MJ/zx547TXo27f284uK4Pe/9+ahb35To37qoSEzgZ8EXjWzM4ETgCV4s05D00REXEmJj/b58ENf26euXbtC8P1+t2+Hn/wE2rdPTj7TRMJ3ACGEUbHHj4DRwELg3BBCeWPSmvRqRKTlKi/3Dt9Fi+Dpp+GrX637PW+8oXb/RmjQWkAhhE+Jj+ZpdJqIZLgQ4Oab4YUX4L//Gy6/vO73fPihT/gaMkTt/g2kTlgRid6kSTB5Mtx6K/zbv9V9ftV2/6uvVrt/AykAiEi0nnkGbr8drrjCA0FdKsf779jh4/3V7t9gCgAiEp3Zs+Haa307xyee8OWdE3nPqlVw2WVq928kBQARicaKFXDJJb6+zwsvQJs2db/nww/hr3/1dv9Ro5o9i+lOAUBEku+jj+D8832i14wZdU/0Am/3nzLFl4JWu3+TUAAQkeQqLPQduoqLE5voBfF2/1271O7fhLQlpIgkz759cPHFiU/0qjRrlrf7X3FF7fsASL0oAIhI8vzwh76ZyzPPJDbRC2D9eu8jGDpU7f5NTE1AIpIcjz0GjzwCP/4xfOMbib1n82Z49FG1+zcTBQARaX5Ll8J3vuP7+d53X2Lv2brVZwUD/OAH0K5d7edLvSkAiEjz2rzZh3v27g1/+hO0TqDleccOL/z37vWZwT0P2jpEmoD6AESk+ZSW+oStwkJf5K1bt7rfU1QE//M/sHOn9xkkMkpIGkQBQESaz623xjt9Tzml7vP37IEHH4QtW7zZ5+ijmz+PGUxNQCLSPJ54Ah5+GG65JbFO33374De/gY8/hhtv9BnC0qwUAESk6eXmwre/7cs0/+IXdZ9fVga//S188AFcfz2cfHLz51EUAESkiW3Z4p2+PXsm1ulbUQFTp8J778FVV0FOTnLyKeoDEJEmVFrqm7kUFMDChdCjR+3nh+C7fy1b5u8bOTI5+RRAAUBEmtJPfgLz58NTT/nM3dqEAH/+s48O+trX4JxzkpNH+Qc1AYlI0/jjH3345s03e1NOXaZPh7lzYfRouOCC5s+fHEQBQEQab9ky+Na3fH2fX/6y7vNffx1efRXOPNM3dNcSD5FQABCRxvniCxg/3tv7p02DrKzaz583zzd1OfVUHx6qwj8y6gMQkYYrKvLO282b4a234PDDaz43BHjxRd8DYPBg3woykS0gpdkoAIhIw6xf7zX/d9+FJ5+sffjm3r3w+OOwfDl85Su+rn+rVsnLq1RLAUBE6u/VV2HiRG++mTEDzjuv5nO3bYPJkyE/HyZMgLPOUrNPitD9l4gkrqICfv5zuPBC35krL6/2wn/jRrj/fp8XcNNNPjNYhX/K0B2AiCRmxw7flGX6dK/9T5lS+968eXne7HPoob6kc+/eycurJEQBQETqtmaNt/evX++rdX7/+zXX5EPwZqGXXoJjj/WNYDp1Sm5+JSEKACJSu+efh2uugQ4dYM4c78StSWmpzwJeuhSGD/cJYYlsACOR0F9GRKpXXg533umreQ4f7oGgT5+az9+50/fv/fBDuPhiGDtW7f0pTgFARA62dasP1Zw929fmf/BBaNOm5vM/+cTX8t+1y8+vax0gSQn1HgVkZq3NbJOZzYsdJ5nZXWb2jplNrnJeQmkikmKWLYNhw+DNN+EPf/B1+msr/Fev9ruEigr48Y9V+LcgDRkGejLwbAhhVAhhFJANnAGcBmwxs3PNbFgiaU1yBSLSNMrLfSP2kSP9+Vtv+eYsNdm711fz/M1v4Igj4I47fGiotBgNaQIaDlxoZmcBq4H3gb+EEIKZvQ6MA3YkmDb7wF9uZjcANwD069evIdckIvX1t7/BddfB4sW+Mudjj9W+rMO6dT77t6DAJ3ZdcglkZycvv9IkGnIH8A5wbgjhNCALaAd8EnutEOgJdEgw7SAhhCkhhJwQQk6PujaTEJHGKSvz5ptTTvEg8PTT8PLLNRf+lbX+X//ah3veeqv3Fajwb5EacgewKoSwN/Y8l3gQAOiIB5XdCaaJSFRWr/Zaf26u1+AnT/amnJocWOsfP772vgFJeQ0phJ82s8Fm1gq4GK/ZnxF7bTCwEchLME1Ekq20FO6+2zt6P/rIl3B+/vmaC/99+/ycA2v9KvxbvIbcAdwN/C9gwHTgHmCBmT0IjI0dHwH3J5AmIsm0fLkvw7xyJVx5pQ/vrK2p9YMPvNa/ZQuMGuV3Cir404aFEBr/S8zaARcAy0IIH9YnrTY5OTkhNze30fkTyXh79/oibv/1X17g//a38PWv13z+vn2+dv+cOdCtm68BNHBg8vIrjWJmeSGEWtbndk0yESyEUAI835A0EWlmS5d6rX/NGvjmN32oZ5cuNZ+vWn/G0ExgkXT1xRfwn//pyzP07g2vvALnn1/7+dOne8Do2hVuuUW1/jSnACCSbvbtg0cegbvu8qUZvvMduPde6Ny5+vN37vTgsGCBb9E4ZozPBVCtP+0pAIikixC8IL/1Vvj7370gf+AB+NKXqj+/pARmzvT1fsrK4IwzvOA/7LDk5lsiowAgkg5Wr/Ymm9mzvdnmlVdg3LjqV+MsLYW5c31z9qIiOPVUuOii2mf+SlpSABBpyQoK4Gc/8925OneGhx6Cb38bsrIOPreiAhYt8pm+27f7ncH48XDkkcnPt6QEBQCRlmjfPnj4YZ/QVVTk++3+x3945+2BQvDx/y++CJs3wzHH+CJvxx2X/HxLSlEAEGlJQvCtFn/8Yx+uef758KtfwaBBB59bUeETvmbM8Bm/vXrBd78LJ5+sjVoEUAAQaRk+/xyefdbH569cCSec4G3455138Lm7dvlSzvPnw7ZtPpHrmmvgy1/2UT4iMQoAIqmqpMRr+0895aN1ysu9w3bKFJ/YdeBeuxs3eudubq6P6hk0yJd7OOkkFfxSLQUAkVRSUeG196eeguee8zH6ffvCT37iG6wf2NRTWgp5eV7wb9wIbdvCmWf6DN7aVvYUQQFAJDWsW+dr8T/9tBfkHTrApZf6GjyjRh1cgy8s9C0bFyyA3bu9ff/KK33z9rZto7gCaYEUAESiUFbmbfnz5/tSzIsXe8fsuef6om3jx3sQOPA9a9fCwoWwYoWnnXKKB4iBA9WxK/WmACCSDHv2wDvvxGvtCxd6zR18PP4vfgETJ0KfPvu/r6zMd+rKy/NCv7gYOnaEsWPhK1+pftinSIIUAESaw65dXquvLPCXLPElmcEL/Kuu8gL8zDMPLvTLy73Qz82NF/rt2nltPycHjj/+4A5gkQbQ/yKRpvDZZz7LduFCL/CXL/eCvFUrGDoUvvc9L/DPOMOHZR6ostDPy/P3Fhd7W35loT9okAp9aXL6HyVSX2VlvvbOokXxY+NGf61NGx9vf8cdXrsfMQI6dar+9+zd64u2rVjhhX5RkRf6gwd7oX/CCSr0pVnpf5dIXbZvh7ffjhf2S5bE2+979YKRI+EHP4DTT4chQyA7u/rfU1HhgWLNGq/tr1/vaW3aeE1/2DAv9Ktbx0ekGSgAiFRVXOw18rw8P955x0fehODNOYMH+6za00/3o1+/mkffhOBr76xd68f773tnsJm/b/Rob9r5p39SoS+RUACQzFVc7EMx8/K8wzUvz2vnFRX++hFHeK38G9/wwv7UU30ETm22b/eCvrLQ377d07t39/cPGuSduAcO8RSJgAKAZI6dO+GZZ3zLw8rCvrzcX+vZ0wv7Sy7xx2HDfBvF2sbW79kDmzbBhg1+bNzoa++AF/DHH+8F/qBBHgBEUowCgKS/8nJ47DH46U99o/PDD/cC/uKL44V9nz61F/bl5fDpp/GCfsMGH/kTgr/eo4c35fTvDwMG1N40JJIiFAAkvc2ZAz/8Iaxa5Z21L7/sTTF11ew//dQL+E8+8QJ/0yZfdwe8dn/00R44+vf3o66mIZEUpAAg6WndOl8z/6WX4KijYNo0X1unasG/d68X8p9+uv9R2YwD3jnbrx989ate0B99tI/jV+1e0oACgKSXbdt8LZ3f/MaHV957ry+otnu3r7JZUBAv6Ldujb+vdWsf0nnccf7Yu7cf3bppKWVJWwoA0vLt2eM1+d/9Dh55xJdh+PKXfWXMjz/2gFCpdWvv8D3mGJ+VW1nY9+ihgl4yjgKApK4QfFOU7du9Zr99e/yo/Lmw0CdVLV7sP/fqBRdd5OvtdO/uBXuPHvHnhx2mgl4kRgFAkq+01Gvp1R07dnjhvm2bN9EUFfn5VY9Wrbwmf8ghvi/uqlW+acovf+lj9tu3j/oKRVoEBQCpv8qa+a5d3qb+xRdeWG/dGi+8d+zwY/duL8RLSnzi1Z49XoiXlx98VFT4Y2mpd9BWTsiqTefOvin6TTd5m7+IJEwBIF1VVHhhW1wcP4qK4o87d/qxbVv8+c6dXqhXFtqV7ysp8WPvXv8FJgZ0AAAFJ0lEQVSdlcsaJyory2vs2dl+tGnji54deqg/tm/vQyvbtfPHTp0SPw49VAumiTRQJN8cM5sKnAC8EkK4J4o8pJQQvOZbVhZv5ti92wvkgoL9a9eFhd7WXVnD3rEj3nxSVOQF9L598THr9XHIIfHCOivLj+xsL6S7d48X0B07euHbubMfhx0GXbr4iJlu3eLt7V27antCkRSW9ABgZpcArUIII8zsMTMbEEJY16Qfsm+fL8JVUbH/UV4eL2wPfK3qOQc2SVR9XvlYVhavYRcV7V+7rqwxVx579sRrzvv2xY+q7dplZfHfW1ZW9zVW1qjbtPGjXTuf4dq2bbyG3a6d166rHh07xo/KQrxLFz86dvT3tG3rR3a2xruLpLEo7gBGAdNiz2cCZwBNGwDmzPENsitVTtc/8PmBaSEcfFRUVJ9e3e+piVm8YK7aBNKpU7ygrix0qxbaHTp47bpbNy+gu3f3o1s3L6yzs+M1dhXUIlJPUQSADsAnseeFwNCqL5rZDcANAP369WvYJ/Tt69P9IT7k75BDvJA08+fVpVeOLmnd2p+3auUFbNXHqke7dvGa9IGPVZ+3aaMCWkRSThQBYDfQLva8I7DfoOwQwhRgCkBOTk49qtlVnHgizJzZiCyKiKS/KGbE5OHNPgCDgY0R5EFEJONFcQfwIrDAzHoD44DhEeRBRCTjJf0OIISwE+8Ifhs4K4SwI9l5EBGRiOYBhBC2ER8JJCIiEdCqWCIiGUoBQEQkQykAiIhkKAUAEZEMZaE+SxokmZkVAB814ld0B75oouy0BJl2vaBrzhS65vo5KoTQo66TUjoANJaZ5YYQcqLOR7Jk2vWCrjlT6Jqbh5qAREQylAKAiEiGSvcAMCXqDCRZpl0v6Jozha65GaR1H4CIiNQs3e8ARESkBgoAacDMWpvZJjObFztOijpP0rTMrKeZLYg972Nm+VX+3nUO95PUZmadzWyGmc00sxfMLDsZ3+m0bALKtE3nzWwoMCGEcFvUeUkGM+sJPB9CONPMsoC/Al2BqSGEx6LNXdMzsy7As8DhIYShsX21e4YQHo04a83CzDoDfwJaAUXABOBR0vg7bWbfBdaFEGaZ2aPAZ0CH5v5Op90dQNVN54FjzGxA1HlKguHAhWa21Mymmlkkq7wmQ6wwfBLfWhTg+0BeCGEkcKmZdYosc82nHC8Ed8Z+Hg78q5ktM7P7ostWs5kIPBBCGAN8DlxBmn+nQwiPhBBmxX7sAZSRhO902gUAqt90Pt29A5wbQjgNyALOjzg/zenAwnAU8b/3m0DaTRYKIew8YN+MGfh1nwqMMLOTI8lYM6mmMPwXMuQ7bWYjgC7ALJLwnU7HAHDgpvM9I8xLsqwKIXwWe54LpF0NqVI1hWEm/r0XhRB2hRDKgeWk6d+7SmH4MRnwNzazrsDDwHUk6TudjgGg1k3n09TTZjbYzFoBFwMro85QEmXi3/t1M+tlZu2BMcC7UWeoqR1QGKb939jMsoHngDtCCB+RpO902v1Dkpmbzt8NPA2sABaHEGZHnJ9kysS/913AXHxb1d+GEN6POD9NqprCMBP+xtcDQ4E7zWwe8B5J+E6n3SggMzsUWAC8QWzTee07nH7MbF4IYZSZHQW8CswGTsf/3uXR5k4aw8y+A9xHvNb7OHAL+k43ubQLAPCPkSKjgTdDCJ9HnR9pXmbWG68hvq6CIT3pO9080jIAiIhI3dKxD0BERBKgACAikqEUAEREMpQCgIhIhlIAEBHJUP8PPwB3CC8AoWcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAD6CAYAAABApefCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XlcVXX+x/HXl0WRRTFEBRWVxHDfcAFTSSs3TNMWmzErm3TKUus3ZU31q5nfVJNtNjU6WppZjjOae7mRhqIXF3DBNRVCUDFZlX39/v4QCBUVLvdyufd+no/Hfdxzv+fecz6n233z9Xs2pbVGCCGEbXOwdAFCCCHMT8JeCCHsgIS9EELYAQl7IYSwAxL2QghhByTshRDCDkjYCyGEHZCwF0IIOyBhL4QQdsDJ0gWUa9asmW7Xrp2lyxBCCKsSExOTqrX2vt376k3Yt2vXjujoaEuXIYQQVkUpdbY675NhHCGEsAMS9kIIYQck7IUQwg7UmzH7qhQVFXHu3Dny8/MtXYrJubi40Lp1a5ydnS1dihDCDtTrsD937hweHh60a9cOpZSlyzEZrTVpaWmcO3eO9u3bW7ocIYQdqNfDOPn5+Xh5edlU0AMopfDy8rLJf7EIIeqneh32gM0FfTlb3S4hRP1U78PeGmVnZ1u6BCGEFdBa891335GYmGj2ddXrMfv6YP/+/bz44os4OTnh4OCAg4MDxcXFODo6AtCyZUuWLVt2zWceeeQRXnnlFUJDQy1QsRDCWpw9e5bw8HDatGmDn5+fWdclPfvb6Nu3L7t27aJjx47MmTOHH3/8kQ4dOrBo0SK2bdt2Q9ADPPvss2zcuNEC1QohrElsbCxKKbp27Wr2dUnPvhq01kRFRfH5559f015UVERubi7/+c9/+Oabbzh16hSdO3eumF/esy8qKuLRRx9lxowZdVm2EKKei42NpUOHDri5uZl9XVYT9itWrCApKcmky2zTpg2PPPLIbd+3ZcsW0tPTGTVqFGfOnOHee+8FYNasWbRt25ZXXnmFadOm0bNnTyIiIkxaoxDCNqWnp5OUlMSECRPqZH0yjHMbOTk5vPrqq6xatYrly5fTp08fAKZOnUqHDh145ZVXLFyhEMIaxcbGAtCjR486WZ/V9Oyr0wM3h7i4OGbNmsVbb72Fm5sbkyZNYuPGjbz00kuMGDHCIjUJIaxfbGwszZs3p0WLFnWyPunZ30b37t158skn+eCDDzh9+jRjx44FIDAwEIDk5GRLlieEsEIFBQX8/PPPdO/evc7WaTU9e0tKTk5m2rRpLF++HEdHRxwdHcnLy6O0tJSnn36ad955h169eqG1ZsCAAdd8trS0FB8fH9atW2eh6oUQ9c3x48cpLi6usyEckLC/rfj4eB5++GG++uqrisOj7r//fqZMmUJJSQkBAQH07NkTgMOHD1uyVCGElYiNjcXV1ZU777yzztYpYX8b/v7+7N69GxcXl4q2CRMm1NkedCGEbSktLeXIkSN07dq14uTMuiBj9tVQOeiFEKI2EhISyMrKqtPxerCCsNdaW7oEs7DV7RJC3FpsbCwODg506dKlTtdbrbBXSrVQSkWWTfsppSKUUtuVUgvVVc5KqQ1Kqd1KqSll77uhraZcXFxIS0uzuWAsv569/ItBCPtz+PBhAgICcHV1rdP13nbMXinVFPgaKD+fdxrwrNb6hFJqE9ANuBeI0Vq/rZTaqJRaCTxzfZvWOqsmxbVu3Zpz586RkpJSo42yBuV3qhJC2I/U1FQuXLhgkfOGqrODtgR4FFgHoLV+vdI8LyAVCAVeLWvbCQTdpO2nygtWSk0FpgJVXvHN2dlZ7uQkhLAZ5WfN1vV4PVRjGEdrfUVrffn6dqXUo8AxrfUFrvb6z5fNSgda3KTt+mUv1FoHaa2DvL29jdwEIYSwDrGxsfj4+GCJvDNqB61Syh/4EzCrrCkbaFQ27V623KrahBDCLuXl5XHq1CmL9OrBiAAuG8NfDkyp1OOPAe4um+4BJNykTQgh7NKxY8coKSmp07NmKzPmpKpXAT/gs7L7qL7F1R24G5VSg4DOwF6uDuFc3yaEEHYpNjYWd3d3i+2HrHbYa61Dy55nA7Ovn6+Uuo+rPfn/1VqXAGeraBNCCLtTWlrK0aNH6datGw4OlhnRNtnlEsp21K64XZsQQtibuLg4cnJyLDaEA7LTVAghzC42NhZHR8drblta1yTshRDCzA4fPsxdd91l0bPmJeyFEMKMLl26xK+//mqxQy7LSdgLIYQZWfKs2cok7IUQwowOHz5Mq1at8PLysmgdEvZCCGEmubm5nDlzxqJH4ZSTsBdCCDM5evQopaWlFh/CAQl7IYQwm9jYWBo3bky7du0sXYqEvRBCmENJSUnFWbNll5axKAl7IYQwg9OnT5OXl1cvxutBwl4IIcwiNjYWJycnAgMDLV0KIGEvhBAmp7UmNjaWwMBAGjZsaOlyAAl7IYQwuYsXL5KSklJvhnBAwl4IIUyu/KzZbt26WbiS30jYCyGEiR0+fBg/Pz+aNm1q6VIqSNgLIYQJZWdnEx8fXy9OpKpMwl4IIUzo6NGjaK0l7IUQwpYdPnwYT09P/Pz8LF3KNSTshRDCRIqLizl27Bjdu3evF2fNViZhL4QQJnLy5EkKCgrq3RAOSNgLIYTJGAwG3N3d6dSpk6VLuUG1wl4p1UIpFVk27ayU2qCU2q2UmlKTNiGEsFU5OTkcPnyY/v374+TkZOlybnDbsFdKNQW+BtzKml4AYrTWA4GHlFIeNWgTQgibtG/fPoqLiwkJCbF0KVWqTs++BHgUuFL2OhRYUTa9EwiqQZsQQtgkg8GAn58frVu3tnQpVbpt2Gutr2itL1dqcgPOl02nAy1q0HYNpdRUpVS0Uio6JSXFuC0QQggLO3fuHImJifW2Vw/G7aDNBhqVTbuXLaO6bdfQWi/UWgdprYO8vb2NKEUIISzPYDDg5OREv379LF3KTRkT9jHA3WXTPYCEGrQJIYRNKS4uZu/evfTs2RM3N7fbf8BCjNll/DWwUSk1COgM7OXqcE112oQQwqYcOXKE7Ozsej2EAzXo2WutQ8uezwL3AbuBe7XWJdVtM3XxQghhabt378bT07NeHltfmVEHg2qtL/DbkTY1ahNCCFtx+fJljh49yogRI3BwqN/nqNbv6oQQoh7bs2cPWmuCg4MtXcptSdgLIYQRtNYYDAY6dOhAixY3HFle70jYCyGEEX755RcuXrxY73fMlpOwF0IIIxgMBho0aECfPn0sXUq1SNgLIUQNFRYWsn//fvr06YOLi4uly6kWCXshhKihgwcPkp+fbzVDOCBhL4QQNWYwGPD29iYgIMDSpVSbhL0QQtRAamoqJ0+eJCQkpN7devBWJOyFEKIGoqKiUEpZxbH1lUnYCyFENWmtiYqKIjAwkKZNm1q6nBqRsBdCiGr6+eefSUtLY+DAgZYupcYk7IUQopoMBgOurq707NnT0qXUmIS9EEJUQ15eHgcOHKBv3744Oztbupwak7AXQohqiI6OpqioyKqOra9Mwl4IIarBYDDg6+tL27ZtLV2KUSTshRDiNpKTk4mPj7e6Y+srk7AXQojbMBgMODg40L9/f0uXYjQJeyGEuIXS0lL27NlDt27daNy4saXLMZqEvRBC3MLRo0e5cuWKVR5bX5mEvRBC3ILBYMDDw4OuXbtaupRakbAXQoibyMrK4vDhw/Tv3x9HR0dLl1MrEvZCCHETUVFRlJaWWv0QDhgR9kqppkqpjUqpaKXUgrK2RUqpKKXUG5Xed0ObEEJYi6KiIn788UcCAwPx9fW1dDm1ZkzP/nFgmdY6CPBQSr0COGqtgwF/pVSAUmr89W0mrFkIIcwuKiqKy5cvM3LkSEuXYhLGhH0a0FUp5Qm0AdoDK8rmbQXuBkKraBNCCKtQWlrK5s2bad++PXfddZelyzEJY8J+F9AWmAGcABoA58vmpQMtALcq2m6glJpaNhwUnZKSYkQpQghhevv37yctLY1Ro0ZZ7Rmz1zMm7N8C/qi1/itwEvgd0KhsnnvZMrOraLuB1nqh1jpIax3k7e1tRClCCGFaWms2bdpEq1at6Natm6XLMRljwr4p0E0p5Qj0B/7Ob8M0PYAEIKaKNiGEqPcOHTpEcnIyI0eOtJlePYCTEZ95D/iKq0M5UcAnQKRSyhcYCQwAdBVtQghRr5X36ps3b06fPn0sXY5J1bhnr7Xep7XuorV211rfp7W+wtUdsnuAe7TWl6tqM2XRQghhDidOnODs2bMMHz4cB4e6OQ2poKCgTtZjkq3RWmdorVdorS/eqk0IIeqzTZs24enpyYABdTMYcezYMbp06cK6devMvi45g1YIIYC4uDhOnTrF/fffj5OTMSPcNbN582ZCQkLIzs7Gx8fH7OuTsBdCCGDjxo24u7tz993mPS1Ia81nn33G6NGjad++Pfv376dfv35mXSdI2AshBElJSRw9epRhw4bRsGFDs62nqKiI6dOnM2PGDMLCwti1axdt2rQx2/oqk7AXQti9zZs34+LiQmhoqNnWkZmZyejRo5k/fz4vv/wyq1evxt3d3Wzru575B6aEEKIe+/XXX4mJiWHEiBG4urqaZR1xcXGEhYURFxfH4sWLeeqpp8yynluRsBdC2LXNmzfj5OTEsGHDzLL8nTt3Mn78eLTWhIeHM2TIELOs53ZkGEcIYbfS09PZs2cPgwYNwsPDw+TL/+qrr7j33nvx9vZm7969Fgt6kLAXQtixrVu3AnD//febdLmlpaXMnj2bKVOmEBoaSlRUFB06dDDpOmpKhnGEEHbpypUr7Nq1i+DgYJo2bWqy5WZnZzNp0iTWrVvHc889x9y5c3F2djbZ8o0lYS+EsEvbtm2juLiYESNGmGyZFy5cYPTo0cTGxvLZZ5/x/PPPm2zZtSVhL4SwO7m5uURERNCnTx+aN29ukmWWn32blpbGDz/8YNI/IqYgYS+EsDsRERHk5+eb7JaDMTExFcsq/yNS38gOWiGEXSkoKODHH3+ke/futG7dutbL2759O/fccw+urq7s2rWrXgY9SNgLIexMZGQkOTk5JunVr1q1ipEjR+Ln58fu3bvp2LGjCSo0Dwl7IYTdKC4uJjw8nI4dO+Lv71+rZX3xxRc88sgjBAUFsXPnTlq1amWiKs1Dwl4IYTeioqLIzMxk1KhRRi9Da827777L1KlTGTFiBOHh4dxxxx0mrNI8JOyFEHYhLy+P9evX4+/vT2BgoFHLKC0t5cUXX+T1119n0qRJrF271mzX0zE1CXshhF3YsGEDWVlZTJw40agbiRcVFTF58mQ+/fRTZs2axddff10vTpaqLjn0Ughh886dO8f27dsZPHgwbdu2rfHnc3JyePjhh9m0aRPvvvsur776qlF/MCxJwl4IYdO01vz73//Gzc2NsWPH1vjz6enphIWFsXfvXhYuXMgzzzxjhirNT8JeCGHT9u7dS1xcHJMnT8bNza1Gn01KSmLkyJGcPn2alStXMn78eDNVaX4S9kIIm5Wbm8t3332Hv78/ISEhNfrsjh07ePjhh8nPz2fTpk0MHTrUTFXWDaN30Cql5imlxpRNL1JKRSml3qg0/4Y2IYSoS+vXryc7O5vHHnus2mPsWmvmzp3LsGHDuOOOO9i7d6/VBz0YGfZKqUFAS631BqXUeMBRax0M+CulAqpqM2HNQghxW0lJSURERDBkyBD8/Pyq9Znc3FwmTZrEiy++yJgxY9i3bx+dOnUyc6V1o8Zhr5RyBr4AEpRSY4FQYEXZ7K3A3Tdpq2pZU5VS0Uqp6JSUlJqWIoQQVdJas3z5ctzd3au9UzY+Pp7g4GCWL1/OO++8w6pVq2jcuLGZK607xvTsJwPHgTlAP2A6cL5sXjrQAnCrou0GWuuFWusgrXWQt7e3EaUIIcSNoqKiiIuLY/z48dU66WnTpk306dOHpKQkNm7cyJ///GccHGzrNCRjtqYXsFBrfRH4FtgJNCqb5162zOwq2oQQwuxyc3NZvXo1d955J8HBwbd8b2lpKX/7298YPXo0bdu2JTo6ut5dh95UjAnhM0D5FYSCgHb8NkzTA0gAYqpoE0IIs1u3bl21dspevnyZ8ePH8+abb/LYY49hMBhqfXG0+syYQy8XAYuVUhMBZ66Oz69XSvkCI4EBgAYir2sTQgizSkxMZMeOHYSGhtKmTZubvu/EiROMGzeOuLg45s6dy4wZM6zujNiaqnHYa62zgIcrtymlQoH7gDla68s3axNCCHOpvFP2gQceuOn7Vq9ezRNPPIGrqyvbtm1jyJAhdVil5ZhkLF1rnaG1XlE2jn/TNiGEMBeDwUB8fDwPPfRQlTtlS0tLeeONN5gwYQJdunQhJibGboIe5AxaIYQNyMnJYfXq1XTo0IH+/fvfMD83N5fJkyezatUqnn76af75z3/SsGFDC1RqORL2Qgirt27dOnJzc6vcKXv+/HnGjh3LgQMH+Pjjj5k1a5bNj89XRcJeCGHVzp49y86dOxk6dOgNNxA/cOAAY8aM4cqVK6xfv56wsDALVWl5cvy7EMJqlV++2MPDgzFjxlwzb/Xq1dx99904OTmxe/duuw56kLAXQlix3bt3k5CQwEMPPUSjRlfP49Ra89577zFhwgR69OjBvn376N69u4UrtTwZxhFCWKVLly6xcuVKAgIC6NevHwAFBQVMnTqVpUuX8thjj7F48WJcXFwsXGn9ID17IYTVKSoqYuHChTg6OjJlyhSUUqSkpDBs2DCWLl3KX//6V5YtWyZBX4n07IUQVmfFihUkJSXx/PPPc8cdd3Ds2DHGjBlDcnIy//3vf3nkkUcsXWK9I2EvhLAq+/btY+fOnQwfPpxu3bqxefNmHn30UVxdXdmxY0fFkI64lgzjCCGsxsWLF/n2228JCAhg7Nix/OMf/2D06NG0b9+effv2SdDfgoS9EMIqFBYWsmDBAho0aMCkSZP4wx/+wMyZMxkzZgy7du265YXPhIS9EMJKLF++nOTkZEaPHs0DDzzAkiVLeOutt1i9ejXu7u6WLq/ekzF7IUS9ZzAYKq43P3HiRHJzc1mzZg3jxo2zdGlWQ8JeCFGvnTt3jmXLlpGens6iRYvw9/fnp59+onPnzpYuzapI2Ash6q38/HzmzZtHZGQkR44cYdSoUSxbtgxPT09Ll2Z1ZMxeCFEvaa35/PPPWbRoEUeOHOH1119n/fr1EvRGkp69EKJemj9/Pm+99RYlJSWsXLmShx56yNIlWTUJeyFEvfP+++/z5z//maZNm7J9+3a5kJkJSNgLIeqNwsJCpk+fzpdffkm7du3YsWMHfn5+li7LJsiYvRCiXjh9+jShoaF8+eWX9OzZky1btkjQm5CEvRDCooqLi/nwww/p3r07sbGxDBs2jA8++ICOHTtaujSbYvQwjlKqBbBZa91LKbUI6Az8oLX+W9n8G9qEEKKyo0ePMmXKFPbv38+QIUPw8/NjwIABDBs2zNKl2Zza9Ow/BBoppcYDjlrrYMBfKRVQVZspihVC2IbCwkL+8pe/0Lt3b3755Rc++OADOnXqRJcuXSquTy9My6iwV0oNBXKAi0AosKJs1lbg7pu0CSEE0dHRBAUF8fbbb/PQQw+xYcMGEhISaNWqFS+88ILccMRMahz2SqkGwJvAq2VNbsD5sul0oMVN2qpa1lSlVLRSKjolJaWmpQghrEheXh6zZ8+mf//+pKWlsX79et59913+/e9/06xZM2bOnImrq6uly7RZxvTsXwXmaa0zy15nA43Kpt3LlllV2w201gu11kFa6yBvb28jShFCWIPIyEh69OjBnDlzmDJlCseOHaNXr158+umneHh4MGvWLDw8PCxdpk0zJuzvBaYrpSKAnsAYfhum6QEkADFVtAkh7ExWVhbPP/88gwcPpqioiPDwcL744gvy8/OZO3cuDRs25KWXXpJLINSBGh+No7UeXD5dFvgPAJFKKV9gJDAA0FW0CSHsyPfff8/06dNJSkpixowZvPPOO7i7u5Oamsonn3yCUooXX3wRLy8vS5dqF2p1nL3WOlRrfYWrO2T3APdorS9X1VbbQoUQ1iExMZEHH3yQMWPG4O7uTmRkJJ9++inu7u5kZmbyySefUFhYyKxZs2jRosrdecIMTHJSldY6Q2u9Qmt98VZtQgjbVVRUxIcffkjnzp3ZsmUL7733HgcPHmTgwIHA1SGdTz75hOzsbGbOnEmrVq0sXLF9kWvjCCFqzWAw8Mc//pEjR44QFhbGZ599Rrt27Srm5+TkMHfuXNLS0pg5c+Y180TdkMslCCGMlpaWxjPPPMPAgQPJzMxkzZo1rF+//powz8/P5x//+AcXL17kueeeIyBAzrG0BAl7IUSNaa1ZsmQJgYGBfPXVV/zpT3/i+PHjjBs37pqzXwsLC/n8889JTExk6tSpcitBC5JhHCFEjRw7doxnn32WyMhIQkJCmD9/fpXXmy8sLORf//oXZ86c4emnn6ZHjx4WqFaUk569EKJacnJyeO211+jZsyfHjh3jyy+/JDIyssqgz8zM5MMPP+T48eM8/vjj9O3b1wIVi8qkZy+EuCWtNWvXrmXWrFkkJiby1FNPMWfOHJo1a1bl+8+ePcu8efPIy8vjueeek7tM1RMS9kKImzp9+jQzZsxg8+bNdOvWjZ07dzJo0KCbvj86OpolS5bQuHFjZs+eLYdX1iMS9kKIG+Tm5vLee+8xZ84cGjZsyNy5c5k+fTpOTlVHhtaa77//nu+//54777yTZ599Vq51U89I2AshKmit2bBhAzNnziQhIYFJkyYxZ84cfHx8bvqZwsJCvv76a6KjowkODmbSpEk3/aMgLEe+ESEEAHFxccycOZMffviBLl26EBERwZAhQ275mczMTObNm0diYiITJkzgvvvukxuP1FMS9kLYuby8PN5//33+/ve/4+zszEcffcQLL7yAs7PzLT939uxZ/vnPf5Kfny87Yq2AhL0QduyHH35gxowZxMfHM3HiRD766CN8fX1v+znZEWt9JOyFsEN79+7lzTffJDw8nMDAQLZt28bQoUNv+znZEWu9JOyFsCOHDh3if//3f9mwYQPNmjXjo48+4vnnn6dBgwa3/WxBQQFLly4lOjqakJAQfv/738uOWCsi35QQduD48eO89dZbfPfdd3h6evLOO+/wwgsvVLtXfvLkSZYuXUp6errsiLVSEvZC2LAzZ87wl7/8hWXLluHm5sabb75Zo9sA5ufns2rVKnbu3Enz5s3505/+RIcOHcxctTAHCXshbNDZs2f5v//7P5YsWUKDBg14+eWXefnll296iYOqHD9+nG+++YaMjAzuu+8+HnjggWoN94j6ScJeCBty4cIF3n33XRYuXIhSiunTp/Paa6/RsmXLai8jLy+PlStXsnv3blq2bMkrr7yCv7+/GasWdUHCXggbEBsby4IFC1i8eDHFxcVMmTKFN954gzZt2tRoOUePHuXbb78lMzOT4cOHM2bMmNseby+sg4S9EFYqNzeXFStWsGDBAvbs2UPDhg353e9+xxtvvFHjnnj5sqKiovDx8eHVV1+VWwfaGAl7IazM0aNHWbhwIUuXLuXy5cvcddddfPzxx0yePBkvL68aLy82NpZvv/2WrKwsRo4cSVhYmBxSaYPkGxXCCpSPoy9YsACDwUCDBg2YMGEC06ZNY/DgwUYdBpmVlcXKlSvZu3cvrVq1Yvr06bRt29YM1Yv6oMZhr5RqAvwHcARygEeB+UBn4Aet9d/K3rfo+jYhRM2cOHGCBQsWsHTpUjIyMujYsSMffvghTzzxRI2OrKksPz+f8PBwwsPDKSoqIiwsjJEjR0pv3sYZ8+3+HvhYax2ulJoPTAQctdbBSqnFSqkAoNv1bVrr06YsXAhbVVhYyJo1a5g3bx47d+7E2dmZ8ePHM23aNEJDQ40+mamoqIiffvqJzZs3k5OTQ+/evXnggQduefliYTtqHPZa63mVXnoDk4C5Za+3AncDvYAV17XdEPZKqanAVAA/P7+aliKETUlKSmLhwoV88cUX/Prrr/j7+/P+++/z5JNP0rx5c6OXW1JSwu7du/nhhx/IzMykS5cujB07VoZs7IzR/25TSgUDTYEE4HxZczrQG3Crou0GWuuFwEKAoKAgbWwtQlir0tJStm3bxrx581i/fj1aa8LCwnjuuee4//77cXBwMHrZWmv279/P+vXrSUlJwd/fn6effpqOHTuacAuEtTAq7JVSdwCfAROAl4BGZbPcAQcgu4o2IUSZjIwMlixZwvz58zl9+jTe3t7Mnj2bqVOn1vqQR601R44cYd26dZw7d47WrVszffp0unXrJtezsWPG7KBtAKwEXtNan1VKxXB1mGYP0AP4GThXRZsQdi8mJoZ58+axfPly8vLyGDhwIG+//TYTJkygYcOGtV7+qVOnWLNmDfHx8TRv3pw//OEPBAUFScgLo3r2T3N1WOZ1pdTrwFfA40opX2AkMADQQOR1bULYpYyMDP773/+yePFi9u/fj5ubG5MnT+bZZ5+lR48etV5+SUkJBw8eJCIigtOnT+Pp6cmkSZMICQnB0dHRBFsgbIHSuvZD5UqppsB9wE6t9cWbtd1KUFCQjo6OrnUtQtQHJSUl/PjjjyxZsoQ1a9ZQUFBA165deeaZZ3jiiSdo0qRJrdeRmZlJZGQkkZGRXL58GS8vL4YOHcqQIUPkEgd2RCkVo7UOut37THJgrdY6g9+OvrlpmxC27tSpUyxZsoSlS5dy/vx57rjjDp555hmeeuopevXqVevhFK01p0+f5qeffuLQoUNorenSpQuPP/44Xbp0qdUOXWHb5CwKIWrp8uXLrFixgiVLlmAwGHBwcGDkyJHMnTuXMWPGmGQsPj8/nz179hAREUFycjJubm7ce++9DB48GG9vbxNshbB1EvZCGKG4uJiIiAiWLFnC6tWrycvLo1OnTsyZM4dJkyaZ7ESlCxcuEBERwZ49eygoKKBt27Y88cQT9O3bV4ZqRI1I2AtRTbm5uWzdupW1a9eyYcMG0tPT8fT05Mknn+TJJ5+kb9++JjnqJSMjgwMHDhATE0NcXBxOTk707duX0NBQuRKlMJqEvRCfs91BAAALYElEQVS3kJ6ezvfff8+aNWvYsmULeXl5eHp6EhYWxrhx4xg9ejQuLi61Xk9GRgYxMTEcOHCAuLg4AFq1asX48eMZOHAg7u7utV6HsG8S9kJcJzExkXXr1rF27Vp27NhBSUkJrVq1YsqUKTz44IMMHjzYJEMoaWlpHDhwgAMHDhAfHw9A69atGTt2LH369KFFixa1XocQ5STshd0rLCwkJiaGbdu2sXbtWmJiYgDo3Lkzs2fP5sEHH6RPnz4mGaJJS0sjJiaGmJgYEhISAGjTpg3jxo2jd+/eEvDCbCTshd3Jzs4mKiqq4hj1vXv3kpeXB0BwcDDvv/8+48aNM8k1ZAoLC4mLi+PEiROcOHGCxMRE4OqF/x588EF69+5dq4ucCVFdEvbC5qWkpLBr166KcD948CAlJSU4ODjQq1cvpk2bxqBBg7j77rtrHbylpaUkJCRw8uRJTpw4QXx8PMXFxTg4OODv78/48ePp06eP0deiF8JYEvbC5qSlpREeHs727duJjIzk5MmTADRs2JABAwbw2muvMWjQIIKDg/Hw8KjVurTWJCcnc+LECU6ePMmpU6fIz88Hrg7P3HPPPQQGBhIQEGCS4+2FMJaEvbB6RUVF7Nmzh61bt7Jlyxaio6PRWtOkSRMGDhzIE088waBBgwgKCqp14BYVFZGUlERCQgLx8fH8/PPPXLlyBQBvb2/69u1LYGAgd911V63/kAhhShL2wirFx8ezZcsWtm7dyvbt27ly5QqOjo7079+ft99+m+HDhxMUFFSrC4GVlpZy4cIFEhISKh7nz5+ntLQUAE9PTwIDAysextzsW4i6ImEvrMKVK1fYsWMHW7ZsYcuWLZw5cwaAtm3bMnHiRIYPH87QoUPx9PQ0avlaa1JSUjh79iy//PILCQkJJCYmUlRUBICrqytt27Zl+PDhtGvXjnbt2hm9LiEsQcJe1Dtaa+Li4jAYDBgMBqKiojhy5Ahaa1xdXbnnnnuYMWMGw4cPJyAgoEaHRJaWlpKamkpycjIXLlwgOTm54lEe7M7Ozvj5+TF48OCKYPf29pZrwgurJmEvLC4vL4+YmJiKcDcYDKSkpADQuHFjgoODmTBhAoMGDSIkJKRa4+7FxcWkpaVVBHr588WLFykuLq54X9OmTfH19aVjx474+vrStm1bfH195eqRwuZI2Is6UVRUxKVLl7h06RK//vorv/76K4cOHcJgMHDw4MGKXnVAQACjRo0iJCSEkJAQOnXqVOW4u9aaK1eukJqaWuUjIyODyvdq8PLywsfHh86dO+Pj41PxMMWlDoSwBhL2ola01pw7d47Y2FjOnz9fEeTloV7+nJGRccNnXVxc6NevH//zP/9DcHAwwcHBFZfrLSgoIDMzk7i4ODIzM8nMzCQ9PZ2UlBTS0tJITU2t+ANRztPTEy8vLwICAvD29sbb2xsfHx9atmwphz0Ku2eSO1WZgtypqv4rLi7m5MmTHDp0iEOHDnHw4EEOHTpEenr6Ne/z9PSkRYsWNG/e/Jpnb29vmjZtiru7O40aNcLDw4OcnJyKMM/MzOTy5ctkZGRUHKtemYuLC97e3jRr1uyGh5eXl1zyV9ilOr1TlbA9WVlZxMbGVgT7oUOHOHLkCAUFBcDVE5S6du3K6NGjueuuu/D396dJkyY4OztTUFBAdnY2WVlZFc9ZWVkkJydXHLZYmaOjI56enjRp0gQfHx86deqEp6dnRVv5tAy5CGE8CXs7UFJSQlJSEqmpqaSkpFQMrZRPp6amkp6eTnp6ekXvujzUAdzc3PD19SUoKAgvLy+aNGmCi4tLxdEpiYmJFdd8Kefq6oq7uzseHh40a9aM9u3b4+Hhgbu7e0V748aN8fT0xN3dXY50EcLMJOzrAa01xcXFFBUVVTxXftysrbCwsMrHpUuXOHPmDPHx8RUnAhUWFla57oYNG+Li4lLx7OnpScuWLXFxcaF58+a0adOG5s2b4+rqWuWjUaNGFdPlIe7m5lark5mEEKZn9WFfvuNOa01paWnFc+Xpqtquf5SUlKC1pqSk5Kbzy58rP27WVlxcTElJyTXP17cVFRVVvN9YRUVFpKenk5qayqVLl7h48SJZWVkAODg40Lp1a4KDg2nfvn3F2LaXlxfe3t54eXnh5uaGi4tLxaNy+Ds5Wf3/HkKIMmb9NSulFgGdgR+01n8zxzr27t3L6tWrzbHoazg6OuLo6IiDg0PF9M3aHBwccHJyokGDBjRq1AgnJyecnJxwdHS85vn6tgYNGuDs7IyTkxPOzs4V0wUFBWRkZJCWlkZKSgqpqan88ssvxMTEcPz48YpDDP39/QkLC6N///7069ePXr16yTi3EAIwY9grpcYDjlrrYKXUYqVUgNb6tKnX06dPH1q3bo1SCgcHBxwcHCqmb/VcHsqVH46OjlXOM2Y8ubS0lIKCgoqhlfLp65/LpzMyMq45AajydE5Ozg3Lb9q0Kf369WPChAkV4S6XzRVC3Iw5e/ahwIqy6a3A3YDJw379+vV89NFHpl7sNcqHd8of5cMwVb0unzZ2aKZRo0b4+vri6+tL7969CQsLw8fHp6LN19cXHx8fPDw8ZKemEKLazBn2bsD5sul0oPf1b1BKTQWmwtU79xjDy8uLzp07G1li9ZUP0ZQPudzudYMGDSoeDRs2vOa5qukmTZrg6+tL48aNJcSFECZnzrDPBhqVTbsDN1xsRGu9EFgIV0+qMmYlY8eOZezYscbWKIQQdsGcV3uK4erQDUAPIMGM6xJCCHEL5uzZrwUilVK+wEhggBnXJYQQ4hbM1rPXWl/h6k7aPcA9WuvL5lqXEEKIWzPrcfZa6wx+OyJHCCGEhcgdGoQQwg5I2AshhB2QsBdCCDsgYS+EEHag3typSimVApw18uPNgFQTlmMNZJvtg2yzfajNNrfVWnvf7k31JuxrQykVXZ3bctkS2Wb7INtsH+pim2UYRwgh7ICEvRBC2AFbCfuFli7AAmSb7YNss30w+zbbxJi9EEKIW7OVnr0QQohbkLC3MkopJ6VUolIqouzRzdI1CdNSSrVQSkWWTbdSSp2r9H3f9hA7Ub8ppZoopTYppbYqpdYopRrUxW/a6odx6uKm5vWJUqo38KjWerala6kLSqkWwHda60FKKWdgNXAHsEhrvdiy1ZmeUqopsBxorrXuXXYv5xZa6/kWLs0slFJNgP8AjkAO8CgwHxv+TSulngNOa63DlVLzgWTAzdy/aavu2Ve+qTngr5QKsHRNdWAAEKaU2qeUWqSUMuuVSy2pLPi+5uotLgFeAGK01gOBh5RSHhYrznxKuBp4V8peDwD+oJQ6oJR613Jlmc3vgY+11vcDF4GJ2PhvWms9T2sdXvbSGyimDn7TVh32VH1Tc1u3H7hXa90PcAZGWbgec7o++EL57fveCdjciTda6yvX3fthE1e3uy8QrJTqbpHCzKSK4JuEnfymlVLBQFMgnDr4TVt72F9/U/MWFqylrsRqrZPLpqMBm+v5lKsi+Ozx+zZorbO01iXAQWz0+64UfEnYwXeslLoD+AyYQh39pq097G97U3Mb9I1SqodSyhEYBxy2dEF1yB6/7y1KKR+llCtwP3DU0gWZ2nXBZ/PfsVKqAbASeE1rfZY6+k1b+39Ie7yp+V+Bb4BDQJTW+kcL11OX7PH7/gvwE1dv7/kvrfXPFq7HpKoIPnv4jp8GegOvK6UigGPUwW/aqo/GUUo1BiKBbZTd1FzudWt7lFIRWutQpVRbYCPwIxDC1e+7xLLVidpQSj0LvMtvvdmvgJeQ37TJWXXYQ8URG/cBO7XWFy1djzAvpZQvV3t+WyQEbJP8ps3D6sNeCCHE7Vn7mL0QQohqkLAXQgg7IGEvhBB2QMJeCCHsgIS9EELYgf8HKTaYcRIOlOQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAD6CAYAAABApefCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl8VNX9//HXSTLZEwgkRDCENSD7FkUEZBMLqFVZxFbU6relttVv++i3VVxRa9X6tfqz+nVBgVq3ihZRxAUVkEVQE5Cw7wECISRAdrLNnN8fWUBIgIQkd5J5P33MI3fOvZn5XMd5e3Luufcaay0iItK8+TldgIiINDyFvYiID1DYi4j4AIW9iIgPUNiLiPgAhb2IiA9Q2IuI+ACFvYiID1DYi4j4gACnC6gUHR1tO3bs6HQZIiJNSnJycpa1NuZs23lN2Hfs2JGkpCSnyxARaVKMMXvPZTsN44iI+ACFvYiID1DYi4j4AK8Zs69OaWkpaWlpFBUVOV1KvQsODiYuLg6Xy+V0KSLiA7w67NPS0oiIiKBjx44YY5wup95Yazly5AhpaWl06tTJ6XJExAd49TBOUVERrVu3blZBD2CMoXXr1s3yLxYR8U5eHfZAswv6Ss11v0TEO3l92DdFZWVlTpcgIk2AtZb3N7/P3uxzmip/XhT25+DVV18lLS2t6vnTTz/N559/Xu22brebiRMn8vDDDzdSdSLSVO3L2ccXu74gPT+9wd9LYX8OBg0axHXXXUdGRgbHjh3j008/ZdSoUWzbtg2Px1O13Y4dO5g6dSq/+tWvSEhI4IYbbmDp0qU/2kZEpFJKRgrGGHq36d3g7+XVs3G8QX5+Pl27dmXlypUEBQVxxx138MILLxAYGMjnn3/Om2++yZQpU5g9ezbx8fF07dqV5557DoDx48fz9ddfM3PmTObMmUPXrl0d3hsR8SYpGSl0jupMeGB4g79Xkwn7eZvmsT9nf72+ZvsW7bmh1w1n3GbLli3ccsstPPjgg+zevZtvv/2Wp59+GpfLhcfjYcmSJYwcOZI777yTyZMnM2TIEAoLC3n00Ud59dVXGT16NHfeeaeCXkR+JLsom305+7i+x/WN8n7nFPbGmFjgfWvt8JPaegPPWmvHGmNcwHygFTDbWjunurb6L7/hXXzxxSxfvpz169czatQobrjhBrp27YqfX/kI2IEDBwgICKCkpISMjAy2bt3KwYMHMcaQnZ1NUFBQ1bYiIpU2ZGwAoG9s30Z5v7OGvTEmCngdCDupzQDPAJWnf94FJFtrHzbGfGKMeQ/41alt1tq8uhZ6th54Q4qJiWHPnj0MHjyYbt260b17d6D8YGxISAjffPMNe/bsYdy4cbz88ssMHTqUiIgIFi5cyIgRI7jvvvscq11EvFNKRgqtQ1vTNrxto7zfuXQ53cBUIPekttuApSc9HwnMq1heDiTW0NYkrVq1ijVr1uByuRg0aBBJSUkkJSWxdOlSAgMDAQgNDeUPf/gDW7du5bHHHiMgIIDAwED+8Y9/0KFDB4f3QES8Sam7lC1ZW+gb27fRzrk5a9hba3OttTmVz40xrYFpwNMnbRYGHKhYPgrE1tD2I8aY6caYJGNMUmZmZt32oIF5PB7+/Oc/c8899wCQnJxMYmIiiYmJjBo1CoDi4mKmTJmCMYZnn32WUaNGsXTpUo4dO0ZaWhp33303O3bscHI3RMSLbDuyjVJ3KX3a9Gm096zLYPKTwL3W2tKT2vKBkIrl8IrXra7tR6y1s6y1idbaxJiYs95oxRG7d++me/fudOvWjdLS0tN69sXFxcyZM4dbbrmFjh07snz5csaNG8fChQsZPXo0mZmZvPzyy0RGRjq9KyLiJdYfWk9QQBDdo7s32nvWZTbOCCCh4k+P/saYx4BkYBjwPtAPWFNDW5PTtWtX5s6dC0BERATLli2rWteyZUtWr16Nx+OpOgi7Z88eJ8oUkSbCWsuGwxvoGdOTAL/GmxBZ63ey1narXDbGLLPWPmCM6QB8YowZDvQEvqV8COfUtmZJs21E5FwdyDvAsePHuKbbNY36vuecUtbakTW1WWv3AmOBVcAV1lp3dW11KdBaW5df83rNdb9E5MxSMlIA6BPbeOP1UI+XS7DWHrTWzjv5YG51bbURHBzMkSNHml0wVl7PPjg42OlSRKSRpWSk0LFlRyKDGvc4nlefQRsXF0daWhreOlPnfFTeqUpEfEducS6p2amNPoQDXh72LpdLd3ISkWZj4+GNWGsb7azZk+nIoohII0nJSKFlcEviIhv/r3qFvYhIIyjzlLE5c3OjnjV7MoW9iEgj2H5kO8VlxY0+C6eSwl5EpBGkZKTg8nfRI7qHI++vsBcRaWDWWjZkbKBHdA9c/q6z/0IDUNiLiDSw9Px0sgqzHBvCAYW9iEiDq7xRSWNe5fJUCnsRkQaWkpFC+xbtiQqJcqwGhb2ISAMqKClg17FdjpxIdTKFvYhIA3LyrNmTKexFRBpQSkYKkUGRdGjh7O1JFfYiIg3E7XGzKXMTfWL7OHLW7MkU9iIiDWTn0Z0cLz3u6CycSgp7EZEGkpKRQoBfAD1jejpdisJeRKShbDi8ge7R3QkKCHK6FIW9iEhDyMjPICM/wyuGcOAcw94YE2uMWVGxHG+MWWaMWWKMmWXKuYwxC40xq4wxt1dsd1qbiIiv2HC4/KxZp6dcVjpr2BtjooDXgbCKpl8Dv7HWjgbaA32Au4Bka+1QYLIxJqKGNhERn5CSkUK7iHa0Dm3tdCnAufXs3cBUIBfAWnu/tXZLxbrWQBYwEphX0bYcSKyhTUSk2SssLWTHkR1e06uHcwh7a22utTbn1HZjzFRgk7X2IOW9/gMVq44CsTW0nfoa040xScaYpOZ4U3ER8U2bMzfjsZ5zCvv9OfsboaI6HqA1xnQG/gT8oaIpHwipWA6veN3q2n7EWjvLWptorU2MiYmpSykiIl4nJSOFsMAwOkV1OuN2e7P3kvB8Av/49h8NXlOtw75iDP8d4PaTevzJwLCK5X5Aag1tIiLNmsd62Hh4I33a9MHPnDliH1j6AMYYrr/o+gavK6AOvzMDiAeerzj9dyblB3A/McYMB3oC31I+hHNqm4hIs7b72G4KSgrOeqOStelreTPlTWYMnUH7Fu0bvK5zDntr7ciKn/cA95y63hgzlvKe/EPWWjewt5o2EZFmLSUjBT/jR6+YXjVuY63lz1/8mdYhrZkxbEaj1FWXnn21Kg7Uzjtbm4hIc5aSkUK31t0IcYXUuM1nOz9jyZ4lPDfuOVoEt2iUunQGrYhIPckqzCI9L/2MQzhuj5u7v7ybLlFduCPxjkarrd569iIivu77A98D0P+C/jVu8/r619l4eCPzJs8j0D+wsUpTz15EpD5Ya/lm/zcktE4gOjS62m0KSwt5cOmDDL5wMJN7Tm7U+tSzFxGpB7uP7eZwwWHGJ4yvcZtnVz/LwbyDvDv53Ua/mYl69iIi9WB12moC/QMZ2HZgtesPFxzmb6v+xnUXXcew+GHVbtOQFPYiIuepxF3C9we+Z2DbgQQHBFe7zaNfP0phaSFPjnmykasrp7AXETlPPxz6gaKyIi5rf1m167cf2c4rya8wfdB0ukd3b+TqyinsRUTO0zf7v6F1aGu6te5W7fp7v7qX4IBgZo6Y2ciVnaCwFxE5D8eOH2Nr1laGxA2p9qDrqn2rmL9lPndfdjex4add/LfRKOxFRM7DmrQ1WGu5NO7S09ZVXhahbXhb/jjkjw5Ud4KmXoqI1NHJc+tjwk6/TPv8LfNZnbaa1655jbDAsGpeofGoZy8iUkeVc+urOzBb4i5hxlcz6BXTi1/0/0XjF3cK9exFROroTHPrZyXPYufRnSz6+SL8/fwdqO7H1LMXEamDM82tzynK4ZGvH2FUx1GM71rzGbWNSWEvIlIHZ5pb/9Sqp8gqzOJ/x/5vo18WoSYKexGROqhpbn1abhrPrHmGm/rcxKB2gxyq7nQKexGRWjrT3PqHlj6Ex3p4bPRjDlVXPYW9iEgt1TS3fkPGBv75wz/570v+m44tOzpTXA0U9iIitXCmufX3LbmPFsEtuG/4fQ5VV7NzCntjTKwxZkXFsssYs9AYs8oYc3tt2kREmrqa5tav3LeSj7d/zIyhM4gKiXKoupqdNeyNMVHA60Dl6V93AcnW2qHAZGNMRC3aRESatOrm1ltrmfHlDNpFtOOuwXc5WF3NzqVn7wamArkVz0cC8yqWlwOJtWj7EWPMdGNMkjEmKTMzs/bVi4g0oprm1n+8/WNW7V/FzBEzCXWFOlhhzc4a9tbaXGttzklNYcCBiuWjQGwt2k597VnW2kRrbWJMzOnXlRAR8SbVza13e9zct+Q+ElolcFv/2xys7szqcoA2HwipWA6veI1zbRMRabKqm1v/9oa32Xh4I4+NfgyXv8vB6s6sLgGcDFTeQLEfkFqLNhGRJqm6ufXFZcU8uPRBBrUdxOSekx2u8MzqciG014FPjDHDgZ7At5QP15xLm4hIk1Q5t35I+yFVba8kv8LenL28es2r+BnvHrw45+qstSMrfu4FxgKrgCuste5zbavv4kVEGkPl3PpurbsRHRoNQF5xHo8tf4zRnUZzRecrHK7w7Op0iWNr7UFOzLSpVZuISFNTObd+fMKJK1g+s/oZMgszeXLMk15zsbMz8e6/O0REvMDqtNUEBQQxqG35hc0yCzJ5evXTTOoxiYsvvNjh6s6Nwl5E5AxOnlsfFBAEwF9X/JXC0kKvu9jZmSjsRUTO4NS59anZqbyU9BK397+di6Ivcri6c6ewFxE5g8q59QmtEgCYuWwmBsPMkTMdrqx2FPYiIjU4dW79xsMbeWP9G9x1yV3ERcY5XV6tKOxFRGqwct/KH82tv3/J/UQGRXLv8Hsdrqz2FPYiItUoKitiyZ4l9LugH9Gh0azat4qPtn3E3UPvplVIK6fLqzWFvYhINVbuW0lhaSHju44vv4TxVzO4IPwCfj/4906XVid1OqlKRKQ5K/OU8cWuL+ge3Z1OUZ1YtH0RK/et5MUJLxIWGHb2F/BC6tmLiJxiTdoasouyGdd1HB7r4d6v7qVLVBd+OfCXTpdWZ+rZi4icxGM9fL7zc+JbxNMjugdvb3ibDYc38PbEt736EsZno569iMhJ1qavrboOTqmnlAeXPkj/C/oztfdUp0s7L+rZi4hUsNby2c7PiA2Ppf8F/Xnp+5fYk72HT2/61OsvYXw2Tbt6EZF6tDlzM/tz9vOTLj+h1F3KEyufYGj7ofyky0+cLu28qWcvIlLh052fEhUSxeC4wcxKnsWBvAP887p/NolLGJ+NevYiIsCuo7vYcWQHYzuPxe1x88TKJ7is/WWM6TTG6dLqhXr2IiLAZzs/IywwjGHxw5j7w1zSctOY/dPZzaJXD+rZi4hwMO8gKRkpjO40GmMMj694nEvjLmVs57FOl1Zvah32xpgoY8wnxpgkY8wrFW2zjTGrjTEPnLTdaW0iIt7os52fERQQxKiOo5i7bi77c/fz8IiHm02vHurWs78ZeMtamwhEGGPuBvyttUOAzsaYBGPMxFPb6rFmEZF6k1WYxfcHvmd4/HBc/i4eX/k4gy8czJVdrnS6tHpVlzH7I0BvY0xLoD2Qw4mbii8GhgEDqmnbcX6liojUv8W7FmOMYWyXsbz+w+vsy9nHy1e93Kx69VC3nv1KoAPw38AWIBA4ULHuKBALhFXTdhpjzPSK4aCkzMzMOpQiIlJ3ucW5rNq3iiFxQwh1hfL4yse55MJLGNd1nNOl1bu6hP1M4A5r7aPAVuDnQEjFuvCK18yvpu001tpZ1tpEa21iTExMHUoREam7L3d/idu6+UnXn/Cv9f8iNTuVmSNmNrtePdQt7KOAPsYYf2Aw8CTlwzQA/YBUILmaNhERr1FYWsjXqV8zqO0gooKj+OuKv5LYLpHxXcc7XVqDqMuY/RPAXMqHclYDzwIrjDHtgPHApYCtpk1ExGt8nfo1RWVFjOs6jjdS3iA1O5Xnxz/fLHv1UIeevbX2O2ttL2ttuLV2rLU2FxgJrAFGWWtzqmurz6JFRM5HqbuUr/Z8Ra82vbgg/AL+uuKvDGo7iKsSrnK6tAZTL2fQWmuPcWL2TY1tIiLeYNX+VeQV5zGu6zje2vAWu4/t5qMbP2q2vXrQ5RJExMe4PW4W71pM56jOdGrZiQlvTWBg24Fc3e1qp0trUAp7EfEp3x/8niOFR7ix9428veFtdh3bxYc3ftise/WgsBcRH1J5c5J2Ee3oEd2DSfMmMeCCAVzT7RqnS2twCnsR8RkpGSmk56Vz+4Db+ffGf7Pz6E4+mPpBs+/Vg8JeRHyE2+NmwdYFRIdGM6DtAG7+4Gb6xfbj2u7XOl1ao1DYi4hPWJa6jIN5B/ntxb9l3qZ57Di6g/k3zPeJXj0o7EXEB+QW5/LRto/o1aYXvWJ6MfX9qfSN7cu1F/lGrx4U9iLiA+ZvmU+pp5SpvaYyb/M8th3ZxvtT3sfP+M79m3xnT0XEJ+06uovV+1cztvNYokOj+cvyv9CnTR+u73G906U1KvXsRaTZ8lgP72x8h6iQKCYkTGDepnlszdrKe1Pe86lePahnLyLN2Mp9K9mfs5/JPSfjZ/yYuWwmfdr0YWKPiU6X1ujUsxeRZqmgpIAFWxfQPbo7g9oO4qWkl9hxdAcf/+xjn+vVg3r2ItJMLdi6gOOlx7mx943kl+TzyNePMKLDCCYkTHC6NEeoZy8izc7e7L2s2LeCMZ3G0C6iHQ8ve5jDBYdZ+LOFPjOv/lTq2YtIs2Kt5d8b/01EYARXd7uaQ/mHePqbp5nScwqXXHiJ0+U5RmEvIs3K6rTV7D62m0k9JxHiCuHRrx+l2F3MX0f/1enSHKWwF5Fmo7C0kPlb5tOlVRcGXziY7Ue2Myt5FtMHTiehdYLT5TlKYS8izcbCbQvJL8nnZ71/hjGG+5fcT4grhIdGPOR0aY5T2ItIs5CWm8bS1KWM6DCC9i3a823at7y/+X3+NORPxIbHOl2e4+oc9saYF40x11QszzbGrDbGPHDS+tPaREQaQuVB2VBXKD/t/lOstdz95d3EhsXyP5f9j9PleYU6hb0xZjhwgbV2oTFmIuBvrR0CdDbGJFTXVo81i4j8SNLBJHYc2cH1F11PWGAYi3YsYvne5cwcMZPwwHCny/MKtQ57Y4wLeBVINcZcC4wE5lWsXgwMq6GtuteaboxJMsYkZWZm1rYUERGKyop4f/P7dGjZgaHxQ3F73Mz4cgYJrRL45cBfOl2e16hLz/4WYDPwFHAJ8DvgQMW6o0AsEFZN22mstbOstYnW2sSYmJg6lCIivu6THZ+QXZTNz3r/DD/jx+vrX2dT5iYeH/M4Ln+X0+V5jbqcQTsAmGWtPWSMeRO4DAipWBdO+f9A8qtpExGpV4fyD/Hl7i8ZGj+UTlGdOF56nIeWPsTgCwczqcckp8vzKnUJ4Z1A54rlRKAjJ4Zp+gGpQHI1bSIi9cZjPby94W2C/IO4/qLya9P/49t/cCDvAE+NfcpnL4tQk7r07GcDc4wxNwIuysfnPzLGtAPGA5cCFlhxSpuISL35dMenbMvaxq39byUiKIIjhUd4YuUTXN3tai7vcLnT5XmdWoe9tTYPmHJymzFmJDAWeMpam1NTm4hIfdh+ZDsLty/k0rhLGRI3BIDHVzxOXkkeT4550uHqvFO9XPXSWnuME7NvamwTETlfucW5vJr8KrFhsfy8z88xxpCancoL37/AL/r9gl5tejldolfSgVMRaTKstcxZN4fjZceZPmg6QQFBADy49EH8jB+PjHrE4Qq9l8JeRJqMT3d+ypbMLdzY+0YujLwQgB8O/cBbKW/x+8G/Jy4yzuEKvZfCXkSahO1HtvPRto8YHDeYoe2HVrXf8+U9RIVEMWPYDAer8366U5WIeL284jxmr51Nm7A2VeP0AF/u/pLFuxbz9yv/Tsvglg5X6d3UsxcRr1Y5Tl9QWsD0QdMJDggGILsom7s+vYsOLTrwu4t/53CV3k89exHxap/t/IzNmZuZ1nda1Zh8ibuEyfMms/PoThZPW1x1oFZqprAXEa+148gOPtz2IRdfeDHD4stPyrfWcsfHd/DVnq/457X/ZFSnUQ5X2TRoGEdEvFJecR6vrX2NmNAYpvWdVjVO//iKx5n7w1weuvwhbu1/q8NVNh0KexHxOtZa5v4wl/yS/B+N07+z4R0eWPoA0/pO4+GRDztbZBOjsBcRr7N412I2Hd7E1N5Tad+iPQAr9q7gFx/+gss7XM5r17ymC53VksJeRLzKzqM7WbB1AYntEhkePxwon2N/3bvX0allJz6Y+oEOyNaBwl5EvEZ+ST6vrX2N6NBobu53M8YYMgsymfDWBPyMH4t+vohWIa2cLrNJUtiLiFfwWA9z180lrzivapy+qKyI6969jrTcND668SO6tOridJlNlqZeiojjrLW8/sPrbDy8kWl9p9G+RXs81sOtC27lm/3f8N6U9xjSfojTZTZp6tmLiKOstby76V3WpK3h2ouuZXiH8nH6+7+6n3mb5vHUFU8xuedkh6ts+hT2IuKoj7Z9xNI9S7myy5WM7zoegFeTX+XJVU/y60G/5k+X/cnhCpsHhb2IOGbxrsV8suMThncYzsQeEzHGsHjXYn6z6DeM6zqOFya8oCmW9URhLyKOWL53Of/Z/B8uvvDiqitZbsjYwOR5k+nVphfvTn6XAD8dVqwvdQ57Y0ysMWZdxfJsY8xqY8wDJ60/rU1EBOC7A9/x9oa36Rvbl9v634af8ePj7R9z5ZtXEhEUwaKfLyIyKNLpMpuV8+nZPw2EGGMmAv7W2iFAZ2NMQnVt9VGsiDR96w+tZ+66uXRr3Y3pg6ZzrOgYN82/iWveuYaY0BgWT1usO041gDr9jWSMGQ0UAIeAkZy4sfhiYBgwoJq2HedTqIg0fVuztjIreRbxLeL5TeJv+M+W/3DXp3eRU5TDIyMfYcawGQT6BzpdZrNU67A3xgQCDwLXAwuAMOBAxeqjwMAa2qp7renAdID4+PjaliIiTcjuY7t58fsXaRPWhim9pnDjf27ko20fcXG7i5lz7Rx6t+ntdInNWl169jOAF6212RVHyfOBkIp14ZQPDVXXdhpr7SxgFkBiYqKtQy0i0gSk5abx/LfPExEYQauQViTOSqTYXczTY5/mD5f+AX8/f6dLbPbqEvZXAKONMb8D+gPxwH5gDdAP2AakUT50c3KbiPigwwWHeW7NcxwvO87a9LUs37ecER1G8NpPX6Nrq65Ol+czah321trLK5eNMcuAnwIrjDHtgPHApYCtpk1EfMzR40f5+zd/Jzk9mTVpawjwC+Dlq17mV4N+hZ/RzO/GdF6TWK21IwGMMSOBscBT1tqcmtpExHccKTzCA0se4MNtH5Ken86EhAm8fNXLVdenl8ZVL2csWGuPcWL2TY1tIuIbdh3dxfSPp/N16teEB4bzxvVvcFOfm3Q2rIN0epqI1Ksvdn3B9I+nk5qdyphOY3hr4lvEhsc6XZbP06CZiNQLay1PrHyC69+9nn05+3hyzJN8cfMXCnovoZ69iJy346XHmfr+VBZuX0hsWCxLbl3CJRde4nRZchKFvYicl82HN3PV21eRmlM+bDP/hvlEBuu6Nt5GwzgiUmez185m0KuDOJB3gL+M+gtf3vKlgt5LqWcvIrVWVFbEbQtu49+b/k1sWCzvTn6XER1HOF2WnIHCXkRqZfPhzVz9ztXsyd7D0PZDmT91Pm3C2jhdlpyFhnFE5JxYa3lt7WsMmDWA9Px0fnfx7/jqlq8U9E2EevYiclZuj5vfLPoNr659lXYR7fjLqL9wW//bdJJUE6KwF5EzOl56nAlvT2BZ6jL6x/bn+QnPMyx+mNNlSS0p7EWkRnuz9zLmX2PYdWwXVyVcxaxrZtEuop3TZUkdaMxeRKr18baPGfDKAPZk7+G+Yffx4Y0fKuibMPXsReRHCkoK+Nuqv/H0N09jscybPI9JPSc5XZacJ4W9iFRZf2g9jy1/jAXbFhAVHMWSW5bQO1a3C2wOFPYiQkFJAe9uepc31r/B8n3L6da6G1/d8pWGbZoRhb2Ij0vJSOFf6//Fin0r+O7Ad4zqOIoPpn5Ai+AWTpcm9UhhL+KjCkoKmLdpHqvTVrMufR1J6Unc1Ocm5lw7h0D/QKfLk3qmsBfxMcVlxXy15ys+3/k5haWFrM9YT1J6EvcMvYfHxzyue8M2U7UOe2NMC+DfgD9QAEwFXgJ6AoustY9VbDf71DYRcU6Zp4yV+1ayaPsicotz6RTVifc2v8e69HU8P/557rzkTqdLlAZUl579TcAz1tovjDEvATcC/tbaIcaYOcaYBKDPqW3W2h31WbiInBtrLd8f/J4Pt35IVmEWXVp1IdQVyt9W/Y3jpcd5/4b3mdhjotNlSgOrddhba1886WkMMA34fxXPFwPDgAGcuNl4ZZvCXqQRWWvZlLmJD7Z8QFpuGnGRcQyNH8rfv/k7KYdTGNNpDM+Pf54eMT2cLlUaQZ3H7I0xQ4AoIBU4UNF8FBgIhFXTVt1rTAemA8THx9e1FBE5xa6ju/hg6wfsOLKD6NBorrvoOt7a8BYPLXuI9pHteW/Ke0zqMUkXMvMhdQp7Y0wr4HlgEvBHIKRiVTjll2DIr6btNNbaWcAsgMTERFuXWkTkhIN5B1mwdQHrD60nMiiSKT2nsO7QOqa+P5VidzH3D7+fe4fdS1hgmNOlSiOrywHaQOA94F5r7V5jTDLlwzRrgH7ANiCtmjYRaQDWWnYd28XSPUtJTk8myD+I6y66Dj/jx+8//z2bMzczvut4nhv3HAmtE5wuVxxSl579f1E+LHO/MeZ+YC5wszGmHTAeuBSwwIpT2kSkHpW4S/juwHcs3bOUtNw0Ql2hXNnlSnq36c1DSx/i3U3v0rFlRz688UOu6XaNhmx8nLH2/EdPjDFRwFhgubX2UE1tZ5KYmGiTkpLOuxaR5i6zIJOv937Nqn2rKCwtJC4yjpEdRzKg7QD+77v/4y/L/4LbupkxdAZ3D72bEFfI2V9UmixjTLK1NvFs29XLSVXW2mOcmH1TY5uI1I21ls2Zm1kx++EtAAALD0lEQVSaupSNhzdiMAxsO5ARHUZwrOgYb214i0nzJpFZmMm13a/l2Z88S6eoTk6XLV5EZ9CKeLHC0kJW71/NstRlHC44TGRQJFclXEVcZBwLty/k6neuZsfRHQT5B/HT7j9l+qDpXNH5CqfLFi+ksBfxMm6Pm21HtpF0MInvD3xPibuELq26cHmHy9matZUnVj7B6rTVGAwjO45kxrAZTOoxSRcukzNS2It4AY/1sC1rG8npyaxNX0tBSQFBAUH0je1LqaeUz3Z+xv1L7qfMU0bvNr352xV/42e9f0b7Fu2dLl2aCIW9iEM81sP2I9tJOpjEuvR15JfkVwV8iCuEJbuX8Kcv/kRucS7tItrxh8F/4OZ+N9M3tq/TpUsTpLAXaUSVAZ98MJl1h9aRV5xXFfA9onuQkpHCS0kvsTZ9LaGuUKb0nMLNfW9mZMeR+Pv5O12+NGEKe5EGVlhayLasbWzO3FwV8IH+gfSN7cugdoPweDzM+WEOf/z8j+SV5NGnTR9eGP8C0/pO0zi81BuFvUg9K3WXsuvYLrZmbWVL5hb25uzFWktQQBC92/QmsV0iXaK6sGDrAn73ye9Yk7aGIP8gpvaeyq8H/ZohcUN0ApTUO4W9yHmy1rI/dz9bMrewJWsLO4/upNRdip/xo3NUZ67udjU9onvQsWVHth/ZzivJr/D6+tfJLsqme+vuPHPlM9za/1ZahbRyelekGVPYi9SSx3o4mHeQ3cd2szVrK1uztlJQUgBAu4h2XN7hcnpE9yChdQLFZcWsTlvNnHVzWJK6hDVpa3D5uZjYYyJ3JN7BiA4j1IuXRqGwFzkDay1Hjh8hNTuVPcf2kJqdyt6cvZS6SwGIComiX2w/Loq+iB4xPcgvyWfF3hXMSp7Fin0rSMlIwWLxN/4MaDuAJ8Y8we0DbqdNWBuH90x8jcJe5CR5xXnszdlbFeyp2ankl+QD4PJ3Ed8inss7XE6nlp3o0KIDWYVZrNq/iheTXmTF3hXsyd4DQKgrlCFxQ3hoxEMMjx/O4LjBhAeGO7lr4uMU9uKTSt2lHMo/xMG8g6Tnp3Mw7yD7c/aTUZBBUVkRxWXFhLhCCHOF4fJzYYyhqKyIbw98y6Idi8gsyCQ9P53somwAYkJjGBY/jDsvuZPh8cPpf0F/XP4uh/dS5ASFvTRr1YV6el46mYWZWGvJKc5hf+5+DuUdYn/ufko9pdW+jr/xp3Voa6JDo4kJjaFnTE9GdBhBYrtEhsUPo1vrbhp7F6+msJcmr6isiKzCrKpHZkEmWYVZHC44XBXqAH7Gj9ahrckvySc9L511h9axP3c/AD2ie3BH4h3ERcYRHRpdFeqVyy2CW+Bnqr3hmkiToLAXr1fqLuVY0TGOHj9abahXjqlXCg4IJiYshrjIOC658BIC/ALYlLmJFftW8OaGN8svS+AfxKhOo7h76N1MSJhA56jODu2dSONQ2ItjPNZDXnEe2UXZZ3wUlhYCYLG4PW48eIgIjCAyKJLosGg6tOxAaEAoQQFBBAUE4fF4KHIXsevoLl5JfoV1h9YB0D6yPdP6TOOqblcxutNoQl2hTu6+SKNS2Eu9sdZSWFpIfkk+ucW5pOensz9nPwfyDpCen87h/MNkFWZx9PhRcopzKCwtxGM9uK0bj/VUPay1WCzWWjzWQ5mnjDJPGW7rrlU9/safy9pfxpNjnmRCwgR6t+mtcXXxWQr7JqbUXUpRWVGNBxJr4va4KSor4njZ8fKfpcernp+8nF+cT35J+eN42XGK3cUUlxWX/3QXU1JWUv7TXUKJu4RSdyklnvLlgpICCksLq16zpnAOd4XTIrgF4YHhBAcEExwQTIgrhFBXKGGuMAL9A3/0cPm5CPQPJCggiJCAkKrtgwOCCQkIqXG5dWhrIoMi6+Nfu0iTp7A/R8VlxeQW55JTnENucW75clHOj9oqn+eW5FLmKTvn17bWUuYpOy14K0P55OXa9m7rg6n8x5Q//Iwffvjh5+eHn/HD3/jjb/yJDIokvkU8MWExXBB2ARdGXkhcRBwdWnagXUQ7YsNjiQ6NJsBP/9mJNLYG/dYZY2YDPYFF1trHGuI9rLWUuEtqDMfqeq/HS49TUFpAQUkBecV55Jee6M1WPgpKCsq3qdjuXHrSLj8XYYFhhLpCcfm5qoYioHy8ubLek5cr1/kZPwL8AnD5ufD3Kw/PQP9AQgJCysPV+GEw+PlV/Kx4fq6MMT/qFYe6Qqt60uGB4VWPiKCIqvHwiKAIwgPDq7YNdYUS5B+koRCRJqjBwt4YMxHwt9YOMcbMMcYkWGt31Pf7/HbRb3k5+eU6/a7B4PJzEeBfHrIuf1dV4Ab4BRAWGEbL4JYE+AWcPqwQEEigX8Vz//JhBn9z9uuNG2MI8Auoevgb//Kffv5Vrx3gF4DL31VVU01tpw53VA51VNfub/wV0iI+rCF79iOBeRXLi4FhwI/C3hgzHZgOEB8fX6c3GRI3hPT89BNh518+IyPQ70TwVf4MDgguX+8fRFhgWFW4V/acKx+Vwehv/Mt70BXLlT3uk39WDWOc0nZqmFcGuuZqi4gTTOVQQr2/cPkQzj+steuNMVcCA621T9a0fWJiok1KSmqQWkREmitjTLK1NvFs2zVkNzMfCKlYDm/g9xIRkTNoyABOpnzoBqAfkNqA7yUiImfQkGP2C4AVxph2wHjg0gZ8LxEROYMG69lba3MpP0i7Bhhlrc1pqPcSEZEza9B59tbaY5yYkSMiIg7RQVMRER+gsBcR8QEKexERH9BgJ1XVljEmE9hbx1+PBrLqsZymQPvsG7TPvuF89rmDtTbmbBt5TdifD2NM0rmcQdacaJ99g/bZNzTGPmsYR0TEByjsRUR8QHMJ+1lOF+AA7bNv0D77hgbf52YxZi8iImfWXHr2IiJyBgr7JsYYE2CM2WeMWVbx6ON0TVK/jDGxxpgVFcsXGmPSTvq8zzrFTrybMaaFMeZTY8xiY8wHxpjAxvhON/lhnMa4z603McYMBKZaa+9xupbGYIyJBd631g43xriA+UArYLa1do6z1dU/Y0wU8A7Qxlo7sOL2nrHW2pccLq1BGGNaAP8G/IECYCrwEs34O22M+S2ww1r7hTHmJSAdCGvo73ST7tmffJ9boLMxJsHpmhrBpcDVxpjvjDGzjTENejE7J1UE3+tAWEXTXUCytXYoMNkYE+FYcQ3HTXng5VY8vxT4pTFmrTHmcefKajA3Ac9Ya68EDgE30sy/09baF621X1Q8jQHKaITvdJMOe6q/z21z9z1whbX2EsAFTHC4noZ0avCN5MTnvRxodifeWGtzT7kc+KeU7/fFwBBjTF9HCmsg1QTfNHzkO22MGQJEAV/QCN/pph72YcCBiuWjQKyDtTSWFGttesVyEtDsej6Vqgk+X/y8v7HW5llr3cA6munnfVLw7ccHPmNjTCvgeeB2Guk73dTD3hfvc/uGMaafMcYfuA5Y73RBjcgXP+/PjTFtjTGhwJXARqcLqm+nBF+z/4yNMYHAe8C91tq9NNJ3uqn/i/TF+9w+CrwB/ACsttZ+6XA9jckXP+9HgKWU3/HtZWvtNofrqVfVBJ8vfMb/BQwE7jfGLAM20Qjf6SY9G8cYEwmsAL6i4j63uv1h82OMWWatHWmM6QB8AnwJXEb55+12tjo5H8aY3wCPc6I3Oxf4I/pO17smHfZQNWNjLLDcWnvI6XqkYVXcwH4Y8LlCoHnSd7phNPmwFxGRs2vqY/YiInIOFPYiIj5AYS8i4gMU9iIiPkBhLyLiA/4/mc0zD5lSoPYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " - 平均潜伏期为：5.40天\n",
      " - 病毒再生基数：2.39\n",
      " - 确诊R2：0.9082\n",
      " - 死亡R2：0.6594\n",
      " - 治愈R2：0.9676\n",
      " - 模型R2：0.8450\n",
      " - 模型总误差：90.5917\n"
     ]
    }
   ],
   "source": [
    "train = wuhan['新增确诊', '新增死亡', '新增出院']\n",
    "train.columns = ['确诊', '死亡', '治愈']\n",
    "train.accumulate(inplace=True)\n",
    "\n",
    "def searchBestParam(seir):\n",
    "    min_loss, max_r2, best_param, likeli_potential = float('inf'),0.0, None, 0\n",
    "    for potential in range(0, 1100, 25):\n",
    "        seir.fit(9000000, potential, 45, 15, 2, train)\n",
    "        loss, r2 = seir.score(9000000, potential, 45, 15, 2, Y=train)\n",
    "        if loss < min_loss and r2 > max_r2:\n",
    "            print('潜在患者：%.4f | R2：%.4f | 误差： %.6f' % (potential, r2, loss))\n",
    "            min_loss, max_r2, best_param, likeli_potential = loss, r2, seir.P, potential\n",
    "    seir.P = best_param\n",
    "    seir.score(9000000, likeli_potential, 45, 15, 2, Y=train, plot=True)\n",
    "    return seir, likeli_potential\n",
    "\n",
    "seir, potentials = searchBestParam(SEIR())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAD6CAYAAABOIFvoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl8VNX9+P/Xyb5DIEDYEVldEEhAcClhU7DggixFkE/9KoI/+QDSCnys/diK1eqHj7ZUrNUPlNaiFUVFtogIKKAoCYoLIFtAwpqQkHWSzMw9vz/OzCRhJzPJTCbvJ4/7mMyZO2fOnWHue+4957yv0lojhBBCAIT4uwFCCCEChwQFIYQQHhIUhBBCeEhQEEII4SFBQQghhIcEBSGEEB4SFIQQQnhIUBBCCOEhQUEIIYRHmL8bcKWSkpJ0hw4d/N0MIYSoVzIzM3O11s0utV69CwodOnQgIyPD380QQoh6RSl1+HLWk9NHQgghPCQoCCGE8JCgIIQQwuOy+hSUUi2Ad7XWtyqlwoH3gCbAIq31Yl+XXelG2O12srOzKSsru9KniguIioqiTZs2hIeH+7spQog6dMmgoJRKBP4BxLqK/hPI1Fr/Tim1Rin1DjDZl2Va66Ir2Yjs7Gzi4+Pp0KEDSqkreao4D601p0+fJjs7m6uuusrfzRFC1KHLOX3kBMYBha77acAy19+fAam1UHZFysrKaNq0qQQEH1FK0bRpUznyEqIBumRQ0FoXaq0LqhTFAkddf+cBLWqhrBql1MNKqQylVEZOTs552ykBwbfk/RSiYapJR3MxEO36O85Vh6/LqtFav6a1TtVapzZrdsm5F0IIEVy0huProPRYrb9UTYJCJnCL6+8bgEO1UFbvrVixgq+++uqi6+Tm5pKfn3/RdQ4dOnROmdPp5HKurX3s2DEqKirOKT948OAlnyuECCBlJyHncyg//5kSX6rJjOZ/AGuUUrcC1wBfYk7/+LKsXjl27Bhz584lJCQEp9PJvHnzmDZtGrfffjvLli1jxowZtG3blrFjxzJnzhxSUlIAePbZZ/n888/Ztm3beet9/vnn6dy5M23atGHr1q2EhoYC8PbbbxMTE8PIkSMBCAkJ4aabbjrn+b/5zW8YOnQo9913X7XylStX0qxZs3PKhRABqjjL3MbV/sCPyz5S0FqnuW4PA0OBrcAQrbXT12W+27y60bJlSxYsWEB+fj5Llizh17/+NQsXLmTBggXs2rWLNm3aAPDcc8/xt7/9DYDjx4+zY8cOhg4dyrJly86p89ChQxw5coRRo0ahtebEiRPk5uaSm5vL3r17adu2red+bm6u53ldu3ZlyJAhDBkyhL179/LnP//Zc79Tp04AzJgxg9WrV1NUdEWDvIQQ/lJ8ECKbQnhCrb+UupzTEIEkNTVVn537aPfu3XTv3t1PLTLWrl3LgQMHmDx5Mn/6059Yu3Yt+fn5FBYW0r59e55//nluvPFGAGw2Gz//+c/54x//SExMDI8//jjdu3dn+vTpuJP9zZs3j9GjR1fbrpSUFOLj46u97unTp/nuu+8896+99lo2b97M3r176devH8uWLaN79+5cf/319O7dmx07dnjam5uby/3333/BbQqE91WIBs9ywu7noXEPaD2ixtUopTK11pcc3VnvEuJd0syZ8M03vq2zZ0/4058uusqiRYvQWrNx40ZeeOEFbDYbycnJdOvWjYyMDNavX8+MGTPo1asXP/30Ew888ABbt24lMzOTd955hxdeeIGQkMoDtwMHDpyzQw4PDyctLa1a2apVq6rdX7t2LcePH2f58uWsWrWKwsJCkpOTATwBAaBfv37MmzfvokFBCBEAbMfAWQFxHevk5YIvKPjB+vXr+fDDD9m1axdz587F6XSyZMkSIiIiiImJoaioiEWLFnHzzTezdu1a3n//fZ544gny8vKYP38+gwYNYtOmTcTExFzytbZs2VLt/tlHeiUlJTzzzDOMHz+ekSNHMnbsWLp06XJOPdHR0dhsNu82XAhR+0pc/QmxHerk5YIvKFziF31taNeuHQ888IDnfmlpKcOGDSMpKYkOHTqwa9cuz2NKKSIiIhg+fDjt27fnnnvu4emnn2bDhg3cdtttREREAGanXVxcTFxcXLXXSk2tfvS3fv16z98rV67khRdeIDY2lpdffplXXnmF7du388tf/hIwAePZZ5/l1ltvJSsri7Zt2/r6rRBC+FpxFkQnQ9ilfzT6QvAFBT/o0qULiYmJTJgwgZMnT5KVlUVMTAylpaUUFxdX+0WemZnJc889R7Nmzfj1r3/NokWL2LFjB//3f/9Hnz59aNHCzN274447WL58Of/xH//Bm2++yeuvv86BAwfOCRKHDx9m4MCBTJgwgYceesgzIklrzRNPPEFaWhotWrRg4sSJ1fIYLVu2jHvuuacO3h0hRI1ZDig9Ak361NlLSpZUH7Esi6VLl3LvvfeyYsUKJk+eTJcuXejZsycdO5pzgXl5eezcuZPx48czevRo3nvvPV555RXS09PZsGGDJyAAjBgxgtWrV3Pq1Cnuu+8+Nm7cyHvvvUeLFi147733+Pjjj7n66quZO3cuy5cv56GHHgLg6NGjLF68mFtuuYX4+HimTZvG0aNHuemmm1i+fDkAe/fu5ejRo/To0aPu3yghxOUrPWICQx0MRXWToOAjt99+O0lJScyZM4du3brRvXt3br75ZpKTk8nMzKRjx47Ex8fzzjvv0Lp1ayZOnMjs2bP5xS9+wbvvvntOf4JSihdffJH169ejtWbGjBls3LiRl156iYSEBJRSLFy4kKSkJB577DG01pSXlzNv3jwKCwtZvnw5TzzxBHFxcTz55JOsWrWK7du3Y7PZ2LBhA/Pnz/fTOyWEuGzFWaBCILZ9nb2kDEn1k/LyciIjI/3djIuqj++rEEHlwCKT4qLTQ15XdblDUuVIwU8CPSAIIfzMWQ6lR+v01BFIUBBCiMBU8hNoS4KCEEIIzPwEFQoxdTt0XIKCHxUXF1/WepIpVYgGqPggxLaFkLq9JK7MU/CR7du389hjjxEWFkZISAghISE4HA5PZtPk5GSWLl1a7Tljx45l9uzZ56SuqEoypQrRADlKwXYCkgfV+UtLUPCRPn36sGXLFh5++GEefvhhUlNTeeihh3jyySc9Se7O9sgjj7BmzZoLBgV3ptQ5c+Zgt9s5ceKEp4N67969jBw50pMhtWrepK5du3pmK9tsNnbt2sXixYs9de7fv58ZM2YwYcIERo4ceU6SPSGEn5UcMrexdX+NdAkKPqS15osvvuDll1+uVm632yktLeXf//43b7zxBnv37uWaa67xPO4OCna7nXHjxjF9+nQA3njjDR599FHAJMMbN25ctUyp7777LnBuptSwsDCWLVt2wUypbhMnTuSDDz6QpHhCBJriLAiNgOhWdf7SQRcUZs6cyTc+zpLas2dP/nQZOZU++ugj8vLyuOOOO9i/fz9DhgzxtKl9+/bMnj2bKVOm0LNnTzZt2nTJ+iRTqhANVHEWxLSHkNA6f+mgCwr+UlJS4kk5cfXVVzN16lQAHn74YYYPH85jjz3ms9eSTKlCBDF7IZTnQpPel163FgRdULicX/S14cCBA8ycOZOnnnqK2NhYJk6cyJo1a5g1axbDhg2rUZ2SKVWIBqj4kLmt4/kJbkEXFPylR48e9OjRg969ezNhwgTuuusu1qxZQ7du3QBz+c2WLVteUZ2SKVWIBqgkC8KiISrZLy8v8xR86Pjx40yZMoW33nqL0NBQQkNDsdlsWJbFgw8+yNdffw2YnXW/fv2qLX379uWuu+6qVp9kShWigdHaNT+hAyjlrzboerWkpKTos+3ateucsrp24MAB3bt3b71z505P2bvvvqv79eun+/Tpo++77z5tWdYV13vkyBG9dOlSbVmWnj59uv7973+vT5486XncbrfrJUuW6EmTJmnLsnRZWZmeMmWKfumll/Tx48er1XXixAk9Z84cXVpaqv/617/qgoKCi752ILyvQjQo5Xla73xK65wvfV41kKEvYx8rWVJ9qKysjKioKH83w2cC5X0VosHI2wHZH0KXRyGqmU+rliypfhBMAUEI4QfFByE8HiKT/NYECQpCCBEItDbzE+Ku8l9/AhIUhBAiMJTngKPEL6ktqpKgIIQQgaA4y9z6aX6CmwQFP3I4HJe1nqTOFqIBKM6CiESIaOzXZsjkNR96/fXXGT58OG3atAFg/vz5XH/99dx+++3nrOt0Ohk1ahS9e/fmd7/73QXrlNTZQjQA2jKZURtdc8lVa5scKfhQSkoKd999NydPniQ/P5+1a9cycOBAfvzxRyzL8qy3b98+xo0bx+TJk+ncuTNjx45l48aN1daBytTZo0aNQmvNiRMnyM3NJTc3l71799K2bVvPfXcKbTCps4cMGcKQIUPYu3cvf/7znz33O3XqBMCMGTNYvXo1RUVFdfPmCCEuzHYCnGV+P3UEyDwFXykuLsayLCIiIoiMjGTq1KnMnDmT7t27s2DBAnJychgzZgyLFi2iXbt25OTk4N6O4cOHU1RUxIYNG1i8eLFnxz1v3jxGjx5dbduqps52Ozt19rXXXsvmzZsvmDrbnSl17dq15ObmXjBLaiC8r0I0CDlb4fjH0P3XEB536fVr4HLnKQTd6aOZ6TP55oSPU2cn9+RPwy6eaG/37t1MmjSJ3/72txw8eJAvv/yS+fPnEx4ejmVZbNiwgbS0NKZNm8bo0aPp378/paWlPP3007z++usMGjSIadOmeQICSOpsIRqM4iwzWa2WAsKVCLqg4C99+vThs88+Y+fOnQwcOJCxY8fSqVMnzxXRjh49SlhYGBUVFZw8eZI9e/Zw7NgxlFKcOXOGyMjIaldPuxhJnS1EELGcUHIYEnv5uyVADYKCUioRWAo0BzK11lOUUouAa4DVWutnXOvVuMwbl/pFX5uaNWtGVlYWN954I126dKFr166A6VSOjo7m888/Jysri2HDhvHqq69y8803Ex8fz8qVKxkwYABPPPFEtfokdbYQDYDtKFj2gOhPgJodKdwPLNVaL1VKvamUmg2Eaq37K6UWK6U6A9fXtExrvc9nW1fHtm7dyrZt27j//vtJSUnxXF3tzJkz3H333QDExMQwc+ZM9uzZwzPPPENYWBgREREsWLDAc/1lN0mdLUQDUHzQzGCO7eDvlgA1G310GrhOKdUYaAtcBSxzPbYOuAVI86KsXrIsi8cff5w5c+YAkJmZSWpqKqmpqQwcOBCA8vJyxowZg1KKl156iYEDB7Jx40by8/PJzs5m9uzZ7NtXGRMldbYQDUBxFkS3NNdQCAA1CQpbgPbAdGA3EAEcdT2WB7QAYr0oO4dS6mGlVIZSKiMnJ6cGTa59Bw8epGvXrnTp0gW73U5KSgoZGRlkZGSwceNGysvLWbx4MZMmTaJDhw589tlnDBs2jJUrVzJo0CBycnJ49dVXSUhI8NSplOLFF19k/fr1aK2ZMWMGGzdu5KWXXiIhIQGlFAsXLiQpKYnHHnsMrTXl5eXMmzePwsJCli9fzhNPPEFcXBxPPvkkq1atYvv27dhsNjZs2MD8+fP9+I4JIbDsUJrt99QWVV3xkFSl1GJgpta6UCk1C/gDMFBrvU0pNQrohtm5v1WTMq31sxd7/UAdkno5LMu67M7kQFBf3lch6q2iA5D1Blw1EeI7XXp9L9Rm6uxE4HqlVChwI/BHKk/73AAcAjK9KAta9SkgCCHqQPFBUCEQ087fLfGoSUfzc8DfMaeQvgBeAjYrpVoBw4F+gPaiTAghgp/WULjHdDCHRvi7NR5X/NNVa/2V1vparXWc1nqo1roQ02G8DXMaqcCbMl9slBBCBLzyHCg/DY0C6xStT85naK3ztdbLtNYnfFFW31XNYXS5mVAv5HwZUu12+zmvcXbepIuRDKlCBICCXWYoakIQBoWGLj8/n0cffZSjR49y4sQJbr/9dvLz8zl+/DjDhw/n0KFDnuDw7LPPeoaFfvnll7z00kv88MMP7N+//5x6n3/+eXbs2MHmzZsZOnQoI0eOpHXr1ixatIi77rqLpk2bcvfdd3P33Xd70lfMnz+fzZs3V6vn7OypK1eu5M0336yNt0IIcbkKdpm+hABIbVGVBAUfSExMZPLkyezbt4/Nmzdjs9nYvXs327Ztw2azsW7dOmw2G5ZlsWLFCm677TYAysrKOHLkCKWlpcyYMaNanVUzpN56663MmTOHrl278vLLLzN16lRWrlxJamoqH3zwAatWrfLMcn7ggQeYNWsWe/fu9dQVEVH9fKVkSBXCz8pPQ9mpgDt1BBIUfCI7O5vk5GQaN27MX/7yF3788Uc2btzIs88+y3fffceiRYv48ssvefvtt2nfvj0xMTGAmYcQERFBnz59ePzxx3E6nZ4633jjDR599FHP/ZiYGL788kvuuecetm3bRt++fWnfvj1Tp06tNgGtadOmrFq1is6dO1+0zRMnTuSDDz7w8TshhLgsBbvMbYCdOoJgTIh3LB3KfNw9EZUMrYZd8OEzZ84wY8YMXn/9dUaNGsWbb75J586dmT9/Pn/4wx/o1asXKSkpDBw4kG7duvHggw+yf/9+cnJyyM/P9yS4u+eee/jVr34FVM+QunTpUl577TW01qSlpTFs2DDuuOMO4uLi6N+/P9nZ2QD86U9/Ys2aNViWxZ133snKlSsJDQ3l22+/ZdiwYViWxZgxY5g8ebJkSBXCnwp2QUwbiGjk75acQ44UfOC6665j5cqVnDp1iu3bt7NixQqcTifz5s1j0aJFJCUlsWDBAp588kkAlixZwpYtW5g0aRJt2rRhy5YtbNmyxRMQzjZ+/Hg2bdpE48aN6du3L61atQLgxIkTJCUledabOXMm69atw263M336dD7++GPS09Pp0aMH6enprFu3jsmTJwOSIVUIv6nIB9vxgDx1BMF4pHCRX/S16ciRI1RUVNCzZ0/mzJnDNddcw5AhQzwZS//4xz+Sm5vLu+++63nOhg0b6NixI/v37692HQWoniG16qS35557jh07dnDw4EGOHDlCu3btzkmdrZS65OxpyZAqhJ8U7Da3AXDpzfORIwUfWbRoES1btiQ8PJydO3eyfv161q9fz1dffVUtM6nb5s2bad26NbNnz+a///u/z3ncnSH1bKGhocTExNC+fXtycnKIjIykf//+1daxLItx48aRlZV1wfYuW7aMESNG1GBLhRBeKdxtEuBFJPq7JeclQcEHzpw5w6FDh+jcuTOhoaH8/ve/9wSFhQsXEhoaCphU1pZlcfDgQWbMmMFzzz1H3759iY+P5/nnn69WZ9UMqe7nuo8Irr32WoqKihg8eDD3338/gwcP9jyvqKiIrVu3MnjwYK666vxJtiRDqhB+Yi+EkiMBe+oIJCj4xLvvvsu4ceMAM7HsySefJC0tjbS0NKZMmUJFRYXnsYqKCmbPns3ChQtp3bo1AH/5y1/45ptv2LZtm6fOqhlSKyoq6NOnD4MGDcLpdPLII4/www8/MHfuXF599VXmzZvnmacQEhLCkiVLmDp1qqeus/sOJEOqEH7iPnWUEJinjqAGWVL9LdCzpNpsNkJDQ8+ZG+BLR48e9QQUgNLSUiIjIz1HJL4SSO+rEEHh4BJwlEKX/6/OX7o2s6SKi4iOjq7VgABUCwhg5jD4OiAIIXzMXmyuxRzAp45AgoIQQtSNwj0mM2oAnzoCCQpCCFE3CndDZBOIau7vllyUBAUhhKhtjlJzLeZG15jMqAFMgkItWbFiBV999dVF18nNzSU/P/+i60jqbCGCQOGPoK2AP3UEEhR84tixY0yaNIlf/vKX3H///Rw6dIhp06bx2muv8etf/5ojR44AMHbsWDIzMz3Pe/bZZxk+fPgF65XU2UIEicLdENHYTFoLdO5JUfVlSUlJ0WfbtWvXOWV1ybIsnZ+fr++8807tcDj0vffeq1esWKFLSkr08OHDtWVZWmut9+/frydPnqy11vrYsWN6wIAB+sknn9Rvv/32OXVmZWXpRx991HP/448/1r/61a/0e++95ym77bbbznlebm6uTk1N1T/++KOnbMCAAeesd9999+nCwsKLbpe/31chgoLDpvW3T2t9NN2vzQAy9GXsY+VIwQeUUnzxxRcMHToUh8NBnz59ePHFF+nfvz+7d+9m4MCBfPnll1x99dW89tpr2Gw2JkyYwAsvvMC4ceP4+9//zqxZs6qdKpLU2UIEicK9oJ0Bm+vobMGXEC89HU74OHV2cjIMu3iivUWLFqG1ZuPGjbzwwgvYbDaSk5Pp1q0bGRkZrF+/nhkzZtCrVy9++uknHnjgAbZu3UpmZibvvPMOL7zwQrUEdpI6W4ggUbgbwuNNqux6IPiCgh+sX7+eDz/8kF27djF37lycTidLliwhIiKCmJgYioqKWLRoETfffDNr167l/fff54knniAvL4/58+czaNAgNm3a5Ln4ztnGjx/Pfffdx8iRI7nmmmto1aoVBw8ePG/q7JkzZzJgwACmT5/O9OnTAUhLSyM9Pb1anZI6W4g64KyAon2Q2DvgRx25BV9QuMQv+trQrl07HnjgAc/90tJShg0bRlJSEh06dGDXrl2ex9xXWxs+fDjt27fnnnvu4emnn2bDhg3cdtttntnQkjpbiCBQtA8sR705dQQy+sgnunTpQmJiIhMmTCAjI4OsrCxiYmIoLS2luLi42i/yzMxMnnvuObKyshgzZgwLFy4kOzubZ555ptrwVEmdLUQQKNwNYbEQ287fLblsEhR8xLIsli5dyr333suKFSuYPHkyXbp0oWfPnnTs2BGAvLw8du7cyfjx4xk9ejTvvfcer7zyCunp6WzYsIEWLVp46pPU2ULUc5YDivZCQjdQ9WdXG3ynj/zk9ttvJykpiTlz5rB48WK6d++Ow+EgPDyczMxMxo0bh91u55133qF169bcc889REdHM2XKFIYMGXJOfVVTZ48ePZqbbrqJ8ePH43Q6mTZtGpZl8corr1BaWsq4ceN45pln6N27tyd1dtUOZEmdLYQfFB8wfQr16NQRSOpsvykvLycyMrJGz5XU2ULUA0feN0cK3X4NIT74bmrtVWe1pM4OcDUNCCCps4UIeJbTpLaI7+qbgGC3w1//ClUGrdQWCQpCCOFrJVngLPPdqaNdu+DUKbjAsHVfkqAghBC+VrALQiMhrqNv6svMhKZNoX1739R3ERIUhBDCl7RlLqgT3wVCfDCW59Qp+OknSEmpkwlwEhT8qLi4+LLWk/TZQtQjhT+a6yc0utY39WVmQmgo9Ozpm/ouQYKCj2zfvp1bbrmFtLQ0Bg0axJAhQ0hLS2Pw4MEMHjyYCRMmnPOcsWPHsmnTpovWK+mzhahn8jIgPAESunhfl90OO3dC9+510p8A1Dx1NvAKMNL19yLgC+DJKo/XuOxiSyCmzq5q8uTJevv27VprrR988EGdlZV1wXU//PBD/fjjj1/wcX+nzw6k91WIeqHstNY7n9L6xCbf1PfNN1o/9ZTWBw96XRW1mTpbKXUrkKy1XqmUGgWEaq37Ax2VUp29KatJewKF1povvvjinJnCdrudgoIC/va3v3HLLbfQvHlz0tLS+N///V+++uor0tLSSEtL4+abb2bBggWe50n6bCHqmbwMM3u5SW/f1JeZCU2aQIcOvqnvMlxxL4hSKhx4HVijlLoLSAOWuR5eB9wC9PKibN+Vtqmq9PR0Tvg4dXZycjLDLiPR3kcffUReXh533HEH+/fv98xUnjlzJu3bt2f27NlMmTKFnj17XvK0EUj6bCHqFcsB+d+YtBbh8d7Xl5NjOpiHDq3TDKs1OVKYBOwCXgD6Ao8CR12P5QEtgFgvys6hlHpYKZWhlMrIycmpQZNrX0lJCXPnzmX58uW89dZbpKSkAPDwww/TqVMnZs+e7VX948ePZ9OmTTRu3Ji+ffvSqlUrgPOmz163bh12u53p06fz8ccfk56eTo8ePUhPT2fdunVMnjwZkPTZQvhUwQ+mg7npJScNX5467mB2q8l4qV7Aa1rrE0qpfwE3AdGux+IwgabYi7JzaK1fA14Dk+biYo27nF/0teHAgQPMnDmTp556itjYWCZOnMiaNWuYNWtWjdsk6bOFqEfyMiCyKcSePxHlFanawRwb6319V6AmRwr7AfeMjFSgA+a0D8ANwCEg04uyeqlHjx788pe/5H/+53/Yt28fd911FwDdunUD4Pjx41dcp6TPFqKesJ2AkiPQJNU3p3p27wabzcxNqGM1OVJYBCxWSv0CCMf0KXyolGoFDAf6ARrYXMOyeuv48eNMmTKFt956i9DQUEJDQ7HZbFiWxYMPPsgf/vAHevXqhdaafv2qb6plWbRs2ZIVK1Z4ykaMGMG4ceMYPnw4zZs3Pyd99ieffOJJnz1t2jTP89zpsxcuXCjps4WoC3kZZqJaoo9O9fihg9njcoYoXWoBEoGxmBFJXpddbAnUIakHDhzQvXv31jt37vSUvfvuu7pfv366T58++r777tOWZV1xvUeOHNFLly7V5eXlOiUlRc+fP187HA49depU/fDDD2uHw6ELCwv18OHDdWZmptZa6+LiYv3Pf/6zWj19+/atdv+vf/2rLigouOhrB8L7KkTAc5Rp/f0ftP7pfd/Ud+qUGYa6ZYtv6nPhMoekSupsHyorKyMqKqpOXqsu0mcHyvsqREA7vR2OroZOD0FMG+/rS0+H7dth1iyf9idcbupsuciOD9VVQIDzp88WQtQxrc2po+iWEN360utfisNhOpi7davzDmY3SXMhhBA1VXoEbCfNMFRfdDDv2uW3DmY3CQpCCFFTeRkmRXaj631Tn7uD+QIDROqCBAUhhKgJRymc+QEa3wChEd7Xl5MDhw/XWYrsC5Gg4EcOh+Oy1pPU2UIEoPyvQTt9N4N5xw6/zGA+mwQFH3r99dc9OYjApLD+6KOPzruu0+lk1KhR/O53v7tonZI6W4gApDXkZUJse4hq7n19Dgd8841fO5jdJCj4UEpKCnfffTcnT54kPz+ftWvXMnDgQH788cdqv+L37dvHuHHjmDx5Mp07d2bs2LFs3LjxnF/6hw4d4siRI4wzkTysAAAgAElEQVQaNYpbb72VOXPm0LVrV15++WWmTp3KypUrSU1N5YMPPmDVqlWkpppfLA888ACzZs1i7969nroiIqof3s6YMYPVq1dTVFRUi++IEEGq+ACU5/nuKCEAOpjdgm5Iavr+dE4U+zhLalwywzpdPH9RcXExnTp1YsuWLURGRjJ16lRefvllIiIi+Oijj/jXv/7FmDFjWLRoEe3ataNTp078+c9/BmD48OF8+umnPPXUUyxevJhOnToBF06dPX/+fLZt28b06dPp2bMnU6dO5fPPP+fbb78FKlNnN29+8V8w7tTZkiVViCt0OgPCYiHBR/N4AqCD2S3ogoK/7N69m0mTJvHb3/6WgwcPenbe4eHhWJbFhg0bSEtLY9q0aYwePZr+/ftTWlrK008/zeuvv86gQYOYNm2aJyCApM4WIiDZC6HoR2h2s2+uwZybazqYhwzxawezW9AFhUv9oq8tffr04bPPPmPnzp0MHDiQsWPH0qlTJ0+W0qNHjxIWFkZFRQUnT55kz549HDt2DKUUZ86cITIy8qIZTcePH899993HyJEjueaaa2jVqhUHDx48b+rsmTNnMmDAAKZPn8706dMBSEtLIz09vVqdkjpbiBrIyzS3TXx0qiczE0JC/N7B7CZ9Cj7UrFkzsrKyiIuLIyUlhb59+5KamkqvXr0YM2YMLVq0wOFwMGzYMNLT02natCnx8fGsXLmSv/3tb0RGRlarz506GyAkJATl+hXx3HPPeY4gLpU6+2IkdbYQV8hyQt4OiOsEEYne11e1gzkuzvv6fECCgg9t3bqVbdu2ER4eTkpKChkZGWRkZLBx40ZPR29MTAwzZ85kz549PPPMM4SFhREREcGCBQto3759tfokdbYQAaboR7AX+a6D2Z0iO9VH9fmABAUfsSyLxx9/nDlz5gCQmZlJamoqqampDBw4EIDy8nLGjBmDUoqXXnqJgQMHsnHjRvLz88nOzmb27Nns21d5NdIRI0awevVqTp06BVA1Ky3XXnstRUVFntTZgwcP9jzPnTp78ODBkjpbCF86nQERjSDeB5eT1xq2bg2YDmY3CQo+cvDgQbp27UqXLl2w2+3nHCmUl5ezePFiJk2aRIcOHfjss88YNmwYK1euZNCgQeTk5PDqq6+SkJDgqVMpxYsvvsj69eupqKigT58+DBo0CKfTySOPPMIPP/zA3LlzefXVV5k3b55nnkJISAhLlixh6tSpnrrO7jvYsGED8+fPr5s3R4hgUH4aig+avgTlg13njz/CiRPws58FRAezm6TOrkOXujzmlZDU2ULUsWNrTZrsbrMg3Mvz/1rD3/4GFRUwbZrpaK5ll5s6O2iOFOpDcPNVQIDzp872ZUCoD++nEHXGXmhGHSXe4H1AANizxxwlDBhQJwHhSgRWa2ooKiqK06dPy47MR7TWnD59uk6vDyFEQDv1mfl133yA93VpDZs2QdOmcL2Psqv6UFDMU2jTpg3Z2dnk5OT4uylBIyoqijZtfHAVKSHqu4p8Mwy1SQpENPa+vt274eRJGDUq4I4SIEiCQnh4+AVH2QghhFdOfmo6lpv/zPu6tIZPPzVHCddd5319tSDwwpQQQgSKslw4sxOa9IHweO/rcx8lBGBfgltgtkoIIQLBqU0QEg7Nb/G+LndfQlJSwB4lgAQFIYQ4P9tJKPgBmt5oMqJ6a9cuOHUqoI8SQIKCEEKc36lNEBIBSTddctVLcvclNGsG117rfX21SIKCEEKcrfQYFOw2ASEs2vv6fvihXhwlgAQFIYQ418mNJhgk9fO+LsuqPEq45hrv66tlEhSEEKKqkp+gaB80uwVCIy+9/qXs2gU5OfXiKAEkKAghRHUnN5pUFk36eF+XZZkRR82bB3xfgpsEBSGEcCvOMkuzWyE0wvv6fvjBXG5zwICAyoR6MRIUhBACzAihkxsgPME3l9p09yU0b14v+hLcJCgIIQRA0X4oOWLSWYT4IAPQ99+bo4S0tHpzlAASFIQQovIoISIREnt5X5/7KKFFC6hn1ySpcVBQSrVQSn3t+nuRUuoLpdSTVR6vcZkQQtSpwj1gOw4tBkCID65L8v33cPp0vepLcPPmSGE+EK2UGgWEaq37Ax2VUp29KfN2g4QQ4opoy4w4ikyCxj64ZrnTWW+PEqCGQUEpNQgoAU4AacAy10PrgFu8LBNCiLpT8AOUnYIWab659vK2beYoYdCgeneUADUICkqpCOC3wFxXUSxw1PV3HtDCy7LzvebDSqkMpVSGXEhHCOEzlhNOboLoFtDIB/MI8vLMvIRu3aBrV+/r84OahMW5wCta6zOu+8WAOzlInKtOb8rOobV+TWudqrVObdasWQ2aLIQQ55G7FcpPQ4vB3v+q1xpWrzazlu+4wzft84OaBIUhwKNKqU1AT2Aklad9bgAOAZlelAkhRO0ryzFXVWt8HSR08b6+776DAwdg8GBISPC+Pj+54sG4WmvPNelcgeFOYLNSqhUwHOgHaC/KhBCidmkLjn5oZi23Gu59faWlkJ4ObdpAaqr39fmRV70qWus0rXUhpsN4GzBQa13gTZk37RFCiMtyeruZqNZyuG8uoLNuHZSVwciR9SLp3cX4YNoeaK3zqRxF5HWZEELUmoozcPITiO8Mja/3vr6sLPjmG7j1VjMMtZ6r3yFNCCGuhNZwdKX5u/UI7zuX7XZYuRKaNIGf/ezS69cDEhSEEA3HmZ1QdACSh0JEI+/r++wzMwx1xAgID/e+vgAgQUEI0TDYi+H4RxDbDpr4oDP45EnYuhV69oSOHb2vL0BIUBBCNAzH1oBlh9Z3en/ayLLMaaOoKLjtNt+0L0BIUBBCBL+C3VCwC5qnQVSS9/VlZEB2NgwbBjEx3tcXQCQoCCGCm8MGx1ZDdEtI6u99fYWF8MkncPXVcL0PRi8FGAkKQojgdmIdOEqhzZ2+SYu9Zo05fTTCB6OXApAEBSFE8Co+CHlfQ7ObzZGCt3bvhj17zNXUEhO9ry8ASVAQQgQnZwVkfwiRTaH5AO/rKy83RwnJydAveDPySFAQQgSnkxvM7OU2d3l/zWWtYe1aKC42qSxCfXAaKkBJUBBCBJ+SI3D6S2ja18xL8FZGhkll8bOfQevW3tcXwCQoCCGCi70YjrwL4Y0gebD39f30kzlK6NzZ9CUEOQkKQojgYTnhp2VmtFH7cRAa6V19RUWwbBk0bgz33huUo43OJkFBCBEctDazlkt+Mv0I3o42cjhMQKiogF/8wsxebgAkKAghgkNeBuRlQvNbzdXUvJWeDkeOwN13Q/Pm3tdXT0hQEELUf8WH4NhaSOgKLQZ5X9+OHaZz+ZZb4JprvK+vHpGgIISo3yryTT9CZFNoO8r78/7Z2bB6tUljMcgHAaaekaAghKi/nBVw+N/mmsvtf+F9x3JxselHSEiA0aPr/aU1a6LhbbEQIjhoDdnvQ9kpaDfGHCl4w+mEd94Bmw3GjYPoaN+0s56RoCCEqJ9OfWpSYiffBvFXe1/funVw+DDceadJZdFASVAQQtQ/Bbvh5CZI7AlJPshDtHMnfPkl3HRTUKbDvhISFIQQ9YvtpDltFNMGWvsgffWxY+YqalddBUOG+KaN9ZgEBSFE/eEohcNvQUikmbHsbaK7ggJ4+22IjW2wHctnk3dACFE/WA5XCotiM9IoPN67+goKYMkSkxL7F78wgUFIUBBC1AOW3RwhFB+C1ndCjJeZSt0BwWaD+++Hlj64AE+Q8PLYSwghapmzwgSEkkPQ9m5I7OFdfWcHhCBPhX2lJCgIIQKXsxwOLYXSI9DmHu8Dwpkz8I9/SEC4CAkKQojA5CyDrH+B7Ri0HQ2Nr/WuPgkIl0WCghAi8DhscOgNM/y03Vho1M27+s6cMaeMyspg0iRo1conzQxGEhSEEIHFUQJZ/4Ty02aUUUJn7+qTgHBFJCgIIQKHvRiy/mEyn7Yf7336CgkIV+yKh6QqpRoppdYqpdYppd5XSkUopRYppb5QSj1ZZb0alwkhGiB7IRz8O1ScgQ4TJCD4SU3mKUwAXtRa3wacAH4BhGqt+wMdlVKdlVKjalrmm80SQtQrFQVwcImZmHbV/RB3lXf15edLQKihKz59pLV+pcrdZsBE4E+u++uAW4BewLIalu270jYJIeqx8jzTh+AsMwEhpo139WVlmRTYliUBoQZqPKNZKdUfSASOAEddxXlACyDWi7LzvdbDSqkMpVRGTk5OTZsshAg0hfvgwOtglcNVk7wLCFrDF1/AP/9pUlY89JAEhBqoUUezUqoJ8BfgXmAW4L4aRRwm0BR7UXYOrfVrwGsAqampuiZtFkIEEG2Z6yGc/BSik82w08gmNa+vogI+/BC+/x66d4e774ZIL6/C1kDVpKM5AngH+C+t9WEgE3PaB+AG4JCXZUKIYOYoNbOUT34KTXrB1Q96FxDy8mDRIvjhBxg8GMaOvbyAUF5uLqwzfTp8/nnNXz/I1ORI4UGgN/AbpdRvgL8D9yulWgHDgX6ABjbXsEwIEaxKs+Gnd8xchDZ3QmIv766HsG8fLF9u6pgwATp1uvj6x4/DmjWwahV8/DGUlJjLbl53nbnAjkBp7f3ZGKVUIjAU+ExrfcLbsotJTU3VGRkZXrdZCFGHtIa8DDiWblJetxsLMV6c79caNm+GjRuhRQtzTeXExHPXsyzYscMEgVWrIDPTlLdtCyNGmGXgwAZxPWalVKbWOvWS6/kiKNQlCQpC1DPOCji2CvK/hfjO0HYUhHmxEy4vh/ffhz17zKUz77wTwsMrH3c4zFHAu++ao4ITJ8yRRP/+lYHguuu8v2JbPXO5QUFmNAshak/5aTj8NpTnQPIgaHardzvj3Fz4979NP8KwYXDjjaY+rc11lv/5T3jzTTh5Eho1Muv8/OcwfDgkJfluu4KYBAUhRO0o2AXZK0CFQoeJ3s1QtizYvh0++cQcFUyaBB06mOsrL10Kb7wB331nHhsxwjx+xx0QEeGzzWkoJCgIIXzLXmj6Dgp2mSuktRsLEY1qXt+JE7ByJRw9ajqSBw2CDRtgyhRYv94EjH794JVXzMijpk19ty0NkAQFIYRvaAtOfwUnN5i/kwdB0k0QUsPdjN0OmzaZCWl2O8TEQHo6TJ1qRg116AC/+Q1MnAhduvhySxo0CQpCCO+VZsPRVWA7AfGdoNUd3s09+PprWLjQnBI6dgyys015kyYwfry5SM4tt0CIXGbe1yQoCCFqzmGDk59AXiaExUH7sZDQ/co7k8vKzASytWvhgw9g/35THhUFAwaYCWZDhsANN0ggqGUSFIQQV05rOPMdHP8InKXQ9EZoMRBCryC1hMNhgsDf/26GjpaXm2DSpo05EviP/zBHA5Kuok5JUBBCXJmyXDi2GoqzTAK71hMhuuXlP3/XLhMI3njDDB1NSoJbb4VmzUwQGDtWho/6kQQFIcTlsRdBzlY4vR1CwqH1CGiScnmnis6cgbfeMsFg+3YIC4OhQ80ksqgo04l8223Qy8u0F8JrEhSEEBdXUQA5WyD/azOqKPEGaDEYwuMu/jyn08wr+PvfzQzk8nIzA/n3vzcprY8eNfMK+vUzS0xM3WyPuCgJCkKI86vIh1NbIP8bcz+xJzS/BSLOk2PIrbDQzB1Yu9b0Exw7ZnISPfQQ3HWXOWLYtQtOnTKnjCQYBBwJCkKI6spPw6nNcOZbQEGT3tDslvNPQNPapKxes8YEgi1bTAdyQoIZLTRunMk59OWXZnRReLjpN+jfX4JBgJKgIIQwyk6ZYFDwPagwaNoXmt1ssppWVVRUeTSwdm3lHIIePeBXvzJ5hm66yRwVfPopLF4swaAekaAgREOmNRQfNGmtC/eYDuSkmyCpf/U+g/x8M3/g7bdNigm7HeLjTWfxU0+ZxHNt2pijhN27TVK6rCyTe0iCQb0iQUGIhqiiwHQc538DFWdMKutmt0JSPwhz7bzPnIEVK2DZMpOK2m6Hq66CGTNM5tGbbqpMOJebCx99ZDKVlpaafoTBgyElRYJBPSNBQYiGwnJC0Y+QtwOKD5ijhLiOkDwEErqZHEWFhfDhv8wRwUcfmUDQvj3MnGnmD6RUGYLqcMC335oL1xw+bGYad+tm1unYUYaW1lMSFIQIdmU5JhCc2WmujxyeAM1/ZkYTRSSaX/Zvv2MCQXq6GTrapg3853+aQNC3b/UdfE6OCQQ7d4LNZvIRDRkCPXtC3CWGqYqAJ0FBiGBUUWD6CAq+h5Ij5poGCV3NSKK4joAy2UeXLDHBoLAQWreGRx4xgeDGGytzDGlthpDu2WP6C44fh9DQyqOCq66So4IgIkFBiGCgNZSdgMIfTTCwuS53HtUMWt5mJpyFxcKRI/CXP5pgsG+fOd8/ZozJMzRgQPVAkJ1tgsCePXD6tClv08bMPL7hBoiN9cumitolQUGI+spyQulhKNhj+goqCswv9pi20HKo6SeIbOo6PfS+CQSffGJ2+AMGwBNPwL33mlFEYGYgHzhggsCePWboaUiIORLo188cGcTHX7RJov6ToCBEfWIvgpJDULgXivaBs8wMI43rCM3TIKGLOSLIy4OPPzejh95+2+zgO3SA//5vc6nKjh1dmU7PmGsXHDxojhzKysycgk6doHt36NwZoqP9vNGiLklQECKQVRRAyWETCEoOm9nGYHb8jbpDfFdz7eOjJyB9M2x5HTZvhu+/N+u5Tw/98pfws59BQQEcOmRyER06ZO6DORXUrZtZrr7aBAbRIElQECJQaA32M1B8yBUIDpv8QwChURDbzmQljW4Lhwtg9VbY/I5JLXH4sFkvPt7MHxg3zuQW6tLFpKc+dAgWLDBHBmCCRYcOcPPN5vRQUpJ0FgtAgoIQ/uMoAdvxKstRc2QAZgJZbHto0gdyFXx7CL7+BjIXw44dZoYxQIsWZuc/a5ZJO920qRkpdOwYZGTAxo1mPXcQuOkmc9usmQQBcV4SFISobVqDvdCMDqoaBOyFletENoGo1lDSBvbkwI79kLnGnO8vdK0XEWFST48ebeYEdO5srktw/LhZPvnErKeU+eV/1VVmmGmHDtC8uQQBcVkkKAjhK1qDowjKc83VySpOm4ljZSfMpDEwO+aIpmBvBNmhsP80fP8TfPeRuUh9SYlZLyrKJJgbO9ac42/d2pwaKigwk8dOnTJL1QDQqpVZkpMr00+IoKC1pqiiiIjQCKLComr1tSQoCHGlnBXmXH95buVScdrcOisq1wsJh4pIOGWHQyWw5yR8cxC+221GB7klJUHXrnDPPdC2rTklFBlpThGVl5sRQQcOmECRlGRGBiUnSwAIMu4df54tz7OcLj3NqeJTHM8/TlFpEfdedy9pXdNqtR0SFIQ4m7PMJImzF5jbijOmA9h931FqjgpKSsypnSIH5NkhpwyOFcFPeXAwBw4chQp7Zb1xcWYoaL9+5nROkyZm1I97wpjn9Z2mD6B9exMEkpJMH0BsrJwCqse01hRXFFNYXkhBeYG5LSsgtySXY/nHOH7mOMWlxdhsNsrKyqgor0CVKUIrQol2/bO3sEPX2m2nBAXRcFh2cBSbsf7uW3sRlORCwSkoPAXFp6C0AGylJq+PzQalZXDGAfnlcNoGp4rhWAHkO6EYsFz1h4aanXeTJtAoGQa4xvjHx5sgEB9vduphYdC4cfUlMdHcJiWZowRRb2itKXeWU1JRQnFFMSV2c5tfms+pwlPkFOaQU5hDfkk+tnIb5eXlVFRUUF5ejrPCSYg9xLPTjyaapJAkWiS0oGVSS5okNqFx48aeJSkpqda3R4KCqJ+0BqsCnKVQVgD5J1079lwoOg3Fp6H0DNgKwHYGys9AebE5FePe2ZeVQYkNSi0oA2yYnXwpUOK67wiD6ASIize/1KOSIKoNdIys3OEnJECjRubXfWioWS8+vnJp1Kj6zl9+8Qc0S1uUOcqw2W2U2kuxOWyUVpRSaCvkTOkZ8oryyC/J50zJGQpsBRSWFVJWUUZFRQV2u92zOB1OInQEka5/UUQRpaJIjkkmKS6JZk2b0TShKQkJCSQmJnp2/PHx8Sg//v+QoCDqnrZcO/Ryc+te3PedZWCVmVtnGdhL4eQRyM6CY4fh1FE4fQJsJWbHXuEAjfnF7nTdlrmWCqAcIBJCYyE0xizhzSA8BsKjzTn5hHBIijA7evcSE2N+1YeEVN6PjTW3MTHVd/zuJSbm3NNBos5Y2qLCWYHdaafcWU65oxyb3UZhqdmh55fkk1+S79mZF9uKKSovIoww0FBSUUJpRSkOh+OcRWuNQhFOOBFEeG4jVASJkYkkRCeQGJdI49jGJMYmkpSQREJCAvHx8Z4lJiaGkAD//xEQQUEptQi4BlittX7G3+1pcLQF2lm5WA7Qjsrbqn+fXWbZQdvNrXvRdnCUg921OMqhosx132Zy9DsscFrgcILdAXbXbakNcvPgdB7knoEc11JmgR1wANFx0KipyesTEgEhkRARDZExEB0LUTGmUzYqypyKiYw0O2r3zt392PkW9w7fvcTGmvIA/yIHKq01DsvhWWx2G0VlRWaHXFZMYVmhOe1SVkxRWRHFZcWUVJRQUu5aykqw2W3YKmyU2csos5dR7iinzF5Ghb3CBADLjtPpNLfaieW0cFpOLMvCsiy0pdFan9UwzI8H19I1oSu92vQiNjyWVpGtiIuMIy4+jvioeOKj40mITqBRTCMSYhKIi4sjJiaGmJgYYmNjiYqK8usve1/ze1BQSo0CQrXW/ZVSi5VSnbXW+3z+Qpbd/BLF/Z9Dm1MQ6Moy933Lcj1umTJtmQXMY9qqXE9b1de3nOa+5TQXIXFUmI5DewU4HeCwm3Kn3ZQ77KbcvThcO1u73TzHXuHaidrNc+yuW4fTVbfD1ON0mB2602le2+lwtcVZvU2Wa6dvuQKBZVUuTtf2Oi1w6soyS1d5zFVuuR/Xrl/o2vxKr3pr4VrPXe5eLHDoynWtsxZ3WXw8OjkZ3ep6uKYNuk0bdNu26Lg4dFgYOjwcIiLQERHo8PBqt5b7sfBwrPBwdGQkVmio+bS12Uk4nc5qtw6nw3PfspXiLCnCOmlhacuUWRYOy4HT6TQ7OYcDp3bicDpwOB04Lad5TDs8ZZZlYbfsWE5z67TMa9kddk+97jqd2onTctWnHdVey2k5PYvDMnW7d4Luxz3tc5e56rNw7SS1a1u0a4epzE7TidPznrjX0VpjYaFx3WqNRlcrr/ZPVf6NovptVWfdvSxVv6Zn/39xl1U9Uqy6nhPzQ8J9616s6uuMfWQsT//q6Ro0Lvj4PSgAacAy19/rgFsAnweFKeNasbss97yPeRPja/J//ErrvOB9XXmjz3pcn72u+29d+bjl2nCtKr9TVb9/ljKLDgUdds6PK7RyfQdVleepKvfPc2tVWceq+tqq+vONIrOofebBbODoRd64C32QytW4S33QwfNj7/Kd9bvofP9xFKryvgKlFcr1ISnXvxAdYm5VSGUZIYS6/oUQQpgKI0yFEapCCVNhhKtwwkPCiQiNICIkgkgVSWRoJNFh0USFRhETFkNMeAzR4dFEh0cTFR5FREQEERERhIeHExoaSlhYGGFhYef9OzQ09KJ/V73fvHnzunm/64FACAqxVH7V84DeZ6+glHoYeBigXbt2NXqR8EZtqCg3Y8PPt5M9Z+d51t/nfG9cOxoNEAJaK9cOTl3g+6Wq7zzd66nKwStVy7VSVN2hmo5JddZ9c+tpg6dO5brvqlcp13qqyg73ItxVo87ZUaqzCqrdV5U36jzrelarcqjtXqfa4ffZr6lU5XquNqkqL6aq/Ktal1JVdl4hqvpz3c9TlfW6X8dTFmJ2bJ4yReVOTylCVSgoCFWhlesqRWhIqFlPmfVCQkIIUWZxP1btfkgIISGm3vDQcELDzE4zNCSUsNAwwkLMzi4sJIzw0HAUisiISLMzDYsgMjzScxsZFum5HxURRURYBNER0eaxcPOcyDDzeHhYuHldVxuD6RSIqLlACArFgDs3bxxwzslbrfVrwGsAqampNfpx/vL/fV3T9gkhRIMRCL1nmZhTRgA3AIf81xQhhGjYAuFI4QNgs1KqFTAc6Ofn9gghRIPl9yMFrXUhprN5GzBQa13g3xYJIUTDFQhHCmit86kcgSSEEMJP/H6kIIQQInBIUBBCCOEhQUEIIYSHBAUhhBAe6pxEUQFOKZUDHK7h05OA8+e6qL+CbZuCbXsg+LYp2LYHgm+bzrc97bXWzS71xHoXFLyhlMrQWqf6ux2+FGzbFGzbA8G3TcG2PRB82+TN9sjpIyGEEB4SFIQQQng0tKDwmr8bUAuCbZuCbXsg+LYp2LYHgm+barw9DapPQQghxMU1tCMFIYQQFyFBoR5SSoUppX5SSm1yLdf7u02iOqVUC6XUZtffrZVS2VU+r0sOCxS1RynVSCm1Vim1Tin1vlIqIhi+T0qpJkqpoUqpJG/qaTBBQSm1SCn1hVLqSX+3xQd6AG9prdNcy3f+bpA3ztqBhiulViqltiql/p+/21YTSqlE4B+YqwoC3Aj8ocrnleO/1l25C+xE6/P3aQLwotb6NuAEMJd6/n1y/Z9bBfQFNiqlmtX0M2oQQUEpNQoI1Vr3BzoqpTr7u01e6geMUEp95frgAyLbbU2cZwf6n0Cm1vpmYLRSKt5vjas5JzAOKHTd7wc8pJTaoZR61n/NqrGzd6K/oB5/n7TWr2itP3bdbQY4qP/fpx7ALK31H4CPgEHU8DNqEEEBc70Gd2rudVRe6a2+2g4M0Vr3BcKBO/zcHm+cvQNNo/Kz+gyodxOKtNaFZ10XZC1mu/oA/ZVSPfzSsBo6z050IkHwfVJK9QcSgY+p598nrfWnWuttSqmfYY4WbqeGn1FDCQqxwFHX33lACz+2xRe+1Vofd/2dAdSrX2pVnWcHGmyfFcDnWusirbUT+Jp6+nlV2YkeoZ5/RkqpJsBfgP9HkHyflFIK8/lqWWwAAAFHSURBVAMrH9DU8DNqKEGhGIh2/R1H/d/uN5RSNyilQoG7gZ3+bpAPBdtnBfCRUqqlUioGuA343t8NulJn7UTr9WeklIoA3gH+S2t9mCD5PmnjUeBb4CZq+BnVqw/TC5lUHj7dABzyX1N84mngDeAb4Aut9Xo/t8eXgu2zAvg9sBFzydlXtdY/+rk9V+Q8O9H6/hk9CPQGfqOU2gT8QD3/Piml5iilJrnuNgb+SA0/owYxeU0plQBsBj4BhgP95FrQgUUptUlrnaaUag+sAdZjfu30c512EX6ilHoEeJbKX9B/B2Yh36eA4RqwsQyIxByJ/hemT+6KP6MGERTA86YNBT7TWp/wd3vEhSmlWmF+5XwkO5vAJN+nwFfTz6jBBAUhhBCX1lD6FIQQQlwGCQpCCCE8JCgIIYTwkKAghBDCQ4KCEEIIj/8ft4jaxsmCe3IAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def forcast(seir, T):\n",
    "    predict = seir.predict(9000000, potentials, 45, 15, 2, T)\n",
    "    plt.plot(train['确诊'], label='确诊(真实)', color='red')\n",
    "    plt.plot(train['死亡'], label='死亡(真实)', color='black')\n",
    "    plt.plot(train['治愈'], label='治愈(真实)', color='green')\n",
    "    # plt.plot(predict['S'], label='易感(预计)', color='blue', alpha=0.5)\n",
    "    plt.plot(predict['E'], label='潜伏(预计)', color='orange', alpha=0.5)\n",
    "    plt.plot(predict['I'], label='确诊(预计)', color='red', alpha=0.5)\n",
    "    plt.plot(predict['D'], label='死亡(预计)', color='black', alpha=0.5)\n",
    "    plt.plot(predict['C'], label='治愈(预计)', color='green', alpha=0.5)\n",
    "    plt.legend()\n",
    "    plt.show()\n",
    "forcast(seir, 30)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
